Estimation of the Complexity of a Finite Mixture Distribution: From Well- to Less Known Methods

Mixture models occur in numerous settings including random and fixed effects models, clustering, deconvolution, empirical Bayes problems and many others. They are often used to model data originating from a heterogeneous population, consisting of several homogeneous subpopulations, and the problem of finding a good estimator for the number of components in the mixture arises naturally. Estimation of the order of a finite mixture model is a hard statistical task, and multiple techniques have been suggested for solving it. We will concentrate on several methods that have not gained much popularity yet deserve the attention of practitioners. These can be categorized into three groups: tools built upon the determinant of the Hankel matrix of moments of the mixing distribution, minimum distance estimators, likelihood ratio tests. We will address theoretical pillars underlying each of the methods, provide some useful modifications for enhancing their performance and present the results of the comparative numerical study that has been conducted under various scenarios. According to the results, none of the methods proves to be a “magic pill”. The results uncover limitations of the techniques and provide practical hints for choosing the best-suited tool under specific conditions.


Aim and Scope
In multiple applications the collected data may be best described by a multimodal probability mass or density function meaning the empirical distribution contains several regions with high probability mass. Mixture models are a powerful mathematical tool that allows for characterizing such heterogeneous populations, which are believed to consist of multiple homogeneous subpopulations. A great multitude of statistical problems can be cast into the mixture model framework: linear inverse and deconvolution problems, random effects models, repeated measures and measurements error models, empirical and hierarchical Bayes, latent class and latent trait models, clustering, robustness and contamination models, hidden mixture structures, random coefficient regression models and many others [43]. A very important class of mixture models is the class of finite mixtures. These models, which assume a finite number of components, have proved to be very useful and flexible enough to model a vast range of random phenomena thus receiving much attention from both theoretical and practical viewpoints.
In some mixture models applications there is no uncertainty about the number of components in the mixture. This is the case where the components correspond to a well-known existing partition of the population. However, on many occasions this situation is far from realistic and practitioners encounter either the lack or complete absence of a priori information about the actual number of mixture components. In such cases, this number has to be inferred from the data along with the parameters of the component densities. Correct identification of mixture complexity may be of primary interest in itself or may be followed by efficient estimation of all parameters. Due to its practical importance the problem of selecting the optimal mixture complexity has been addressed in numerous statistical publications, and we will point out many seminal works as we proceed.
The objectives of the present survey are to • provide the theoretical background of the reviewed methods for estimating the complexity of a finite mixture; • assess the performance of these methods under various scenarios; • suggest modifications that enhance the performance of some of the methods in particular settings; • identify universal methods that provide stable and accurate results throughout most of the scenarios for different distribution families or single out scenarios under which certain approaches may be preferred to others.
The number of methods devoted to estimating the true number of components in a mixture is undoubtedly too large to be thoroughly described in a single survey. Thus, we restrict attention to a selected subset of approaches that have the merit of being applicable in very general settings; i.e., for wide classes of finite mixtures of distributions. One of the main goals of this survey is to uncover the extent to which each of the methods is successful in consistently estimating the true complexity for various sample sizes. The estimation techniques reviewed below can be split into three main categories: 1. methods built upon the determinants of the Hankel matrix of moments of the mixing distribution; 2. methods based on penalized minimum distance between the unknown probability density and a consistent estimator thereof. The distances considered in this survey are the Hellinger as well as the L 2 -distance; 3. likelihood ratio test (LRT) -based techniques.
Some of the key criteria we based our choice of the techniques upon were: a) a cohesive mathematical theory behind the method, including the asymptotically consistency; b) infrequent reference in the literature as well as relatively rare usage in practice despite of the coherent theoretical base; c) feasibility of implementation using any programming language, e.g. Python, R, Julia, Matlab, etc.
Pseudocodes for the algorithms discussed in this work are given in Appendix D in the supplementary materials. For completeness, other interesting methods for estimating the complexity of a finite mixture are mentioned in Sect. 8.
Although not strictly a part of a survey, the performance enhancement such as the one we bring to some of the methods through specific modifications is almost inevitable. In fact, the original version of some of the approaches reviewed here cannot be of real practical value without any further adjustment. These modifications, which will be described in separate subsections, include resorting to some judicious scaling in the case of the Hankel-based-methods or using bootstrap instead of penalization for the approaches based on minimum distance estimation. In Sect. 6, we report the results of an extensive numerical study which we carried out for different mixture distributions and various number of components with the goal of comparing the performances of the techniques reviewed in this survey.
Several examples involving the estimation of mixture complexity for real data sets using the discussed methods are presented. The data sets were taken from various fields such as geology, insurance and lexicography.

Organization of the Paper
The paper is organized as follows: • Section 2 provides some basic background on mixture models, mentions major works on mixture model estimation techniques and gives a brief overview of these approaches. • Section 3 outlines the theoretical foundation of the original method based on the determinants of the Hankel matrix of moments of the mixing distribution as proposed in [21]. In the same section, two modifications of this approach allowing for obtaining improved results, are presented. The section also gives a concise description of a neural network extension of the Hankel matrix approach, proposing possible working configurations of a multilayer perceptron for the mixtures of Gaussian, Poisson and geometric densities. • Section 4 describes methods based on minimum distance estimation. The section relies to a very large extent on the works [74], [69] and [20]. We re-examine the estimation techniques that use the Hellinger and the L 2 distances when combined with two different penalties. In the same section, motivated by the idea of enhancing the original method, we propose a modification based on a bootstrap procedure instead of penalization. • Section 5 presents the estimation approach based on the LRT combined with a bootstrap procedure as described in [38]. • Section 6 comprises the results of a comparative numerical study where all of the above mentioned techniques are tested on simulated data under various scenarios. Furthermore, the same section contains a discussion of the settings in which certain methods can be favored as they seem to outperform their counterparts. • Section 7 encompasses several real data sets that were analysed using the studied approaches and compares the obtained results. • Section 8 mentions a number of papers where other techniques for mixture complexity estimation, not addressed in the present survey, are considered. • Section 9 summarizes the findings and outlines the limitations of all methods reviewed in this survey. • Appendices A, B, C, D presented as supplementary materials to this paper include additional examples, tables with detailed simulation results, proofs of the theoretical results that are relevant for the methods described in the manuscript and pseudocodes clarifying and simplifying the implementation of the discussed techniques.

Notation and Basic Definitions
We start with defining the terminology that will be used throughout the survey. In the sequel, a real vector of dimension r will be denoted v r and its components by v 1 , . . . , v r . A class of real vectors of dimension r will also bear the subscript r in its notation. When manipulating several vectors of dimension r we will index them as v r ,1 , v r ,2 , . . .. A random sample of i.i.d. random variables are going to be denoted (X 1 , . . . , X n ). Also, a sequence of random variables (for example converging weakly) will be denoted for example by Y (n) . A class of densities which depend on some vector of parameters of dimension r will not necessarily bear the subscript r . Suppose that some population of interest, represented abstractly by a random variable X , consists of a finite number m ∈ N of subpopulations. Each subpopulation is generated by some random process that can be modeled by an individual or component distribution, e.g. normal, exponential, Poisson, geometric, etc. We will assume that each of the component distributions admits a density with respect to some common dominating measure μ. Furthermore, the component density is assumed to be parametrized through some unknown vector φ d ∈ ⊆ R d , d ≥ 1. To keep the manuscript to a reasonable length, we confine our attention in this survey to the onedimensional case; i.e. the random variable X ∈ R. In the sequel, the dominating measure μ is either the Lebesgue measure in case the distribution of the components is absolutely continuous, or the counting measure in case this distribution is discrete.
In the latter case, all the examples considered here treat distributions that are supported on the set of non-negative integers. Let F = f φ d : φ d ∈ be the family of densities which the components belong to. If X is the support of X , then the distribution of X is said to have a m-component mixture distribution with density Above, G is a discrete distribution defined on with at most m jump points at φ d,1 , . . . , φ d,m , and the integral representation in (2.1) is given here only to draw attention that finite mixtures are part of a much bigger family of mixtures where G can be any distribution function, known often under the name of "mixing distribution". In the sequel, we will refer to either the probability density or probability mass function defined in (2.1) as the mixed density and to π j , j = 1, . . . , m as the mixing probabilities.
We define the family of m-component mixture densities as the set where f θ pm is given by (2.1). Suppose we observe n random variables X 1 , . . . , X n ∈ X which are i.i.d. according to an unknown density ∼ f 0 ∈ m≥1 F m . What is the value of m that can be assigned to this density based on the observed data? It is clear that such a value needs to target the most parsimonious representation of the mixture. Estimation of the true complexity of f 0 cannot be presented without touching upon this point, which is discussed in the next section.

Identifiability and Complexity of Mixture Models
We will now touch upon the identifiability issues arising within the mixture distributions framework, which is a crucial point when the aim is to estimate the true complexity. Identifiability of general mixtures of some additively closed family of distributions was proved in the pioneer work of Teicher [65], who recognized the importance of settling the issue of identifiability before launching into estimation of the mixing distribution. Several articles have been devoted to proving identifiability of finite mixtures of some particular classes of distributions such as finite mixtures of normal or gamma distributions; see e.g. [66]. For identifiability results in other classes or review papers on the subject we can refer to [19,35,36,43,49,68].
The identifiability of a mixture is defined as follows: a finite mixture with respect to the family F is said to be identifiable if for any m ≥ 1 and any two elements f θ m and f θ m in F m satisfying the equality Different techniques have been developed to show identifiability. One of the most important results is the one shown in [76], which says that the characterizing condition of identifiability is linear independence of the family F. Other characterizations or sufficient conditions could be built upon this result by resorting for example to using some additional properties of the elements of F or computing Fourier/ Laplace transforms (see [5,36,39]).
When identifiability holds, it is natural to think of the most economic representation of the finite mixture under study. Indeed, we have the inclusions for all m ≥ 1, and hence we can introduce the following definition: the index of economical representation for some finite mixture density f ∈ m≥1 F m is defined as This index is exactly what is called the complexity (or order). Note that this number has to be unique, an immediate consequence of identifiability. Also, from a practical point of view, m( f ) corresponds to the number of all the components that are actually part of the total population: all the mixing probabilities π j , j ∈ {1, . . . , m( f )} should satisfy π j > 0 by the very definition of m( f ). The term identifiability is used here with some abuse as the components of θ p m( f ) are unique up to some permutation (whereas the mixed density is invariant under the m! permutations of the component labels). One can of course require for example that the mixing probabilities are labeled so that π 1 < · · · < π m in case they are all different. We will be following this convention when reporting the simulation results in Sect. 6.
The discussion above lays the ground for this survey. In the sequel, we shall assume that identifiability assumption holds. Also, the notation m( f 0 ) = m 0 will be used, where f 0 is the unknown density in F m 0 from which we observe a random sample. The true complexity or order, m 0 , as well as the true parameter vector θ 0 := θ p m 0 ∈ p m 0 will be assumed to be unknown. The main goal of the methods reviewed further is to consistently estimate m 0 . An estimation procedure can be (but does not necessarily have to be) accompanied by the estimation of θ 0 .

Popular Approaches to Mixture Model Estimation
Mixture model estimation has a long history. The early mixture model estimation techniques date back to the end of the 19-th century, when S. Newcomb [52] suggested an iterative reweighting scheme to compute the Maximum Likelihood (ML) estimator of the common mean of a mixture of a known proportions of a finite number of univariate normal populations with known variances. This scheme is regarded by many as a precursor of the well-known Expectation-Maximization (EM) algorithm.
A few years later K. Pearson [56] described an analytical and a graphical solutions to estimating the first five moments of an asymmetrical empirical distribution, which he was aiming to break up into two univariate normal curves. The graphical solutions for mixture model estimation stayed in the focus of attention until the second half of the 20-th century ( [14,34,58]).
Between 1912 and 1922 R. Fisher [29] attempted to popularize the ML approach to fitting the mixtures. The evolution of the ML approach is considered in detail in [3]. In particular, Fisher made an analysis of the extensions of the method of moments to the likelihood equations as a way of increasing the quality of the estimates, which later caused a dispute with Pearson ( [28,57]). Around the 1950s C. R. Rao [59] used Fisher's scoring method to estimate the parameters of a mixture of two Gaussian distributions with common variance, and soon after the ML estimation for identifying the number of components as well as for parameter estimation in finite mixture models was addressed in numerous publications, such as [22,72,73].
These days the most well-studied and widely-used approach to computing ML estimates for finite mixture models as defined in (2.1), is the EM algorithm, elaborately described in [23], the seminal work that greatly exhilarated the efficient usage of mixture models. The EM algorithm is implemented by assuming that there are latent variables that link every observation to one of the components, which, together with the observed data, yield complete data.
We will summarize the main idea behind this algorithm. To that end consider two sample spaces within the mixture model framework: 1. the sample space of the incomplete observations, where only the realizations of the random variable X are observed, but no information on the mixing distribution G(φ d ) is available; 2. the sample space of the complete observations, the estimation of which can be performed explicitly.
For the sake of simplicity consider the one-dimensional case, ⊆ R. In this case we denote φ d simply by φ. The extension to the multidimensional case is possible but complicates the derivations. Let x = (x 1 , . . . , x n ) be the observed realizations of the random variable X , and let z = (z 1 , . . . , z n ) denote the realizations of the corresponding unobserved (or latent) random vector Z indicating that the observation x i , i = 1, . . . , n comes from the jth component, j = 1, . . . , m. In other words, z i , i = 1, . . . , n are realizations of a multinomial distribution with probabilities π 1 , . . . , π m , and we have that The pairs y i = (x i , z i ), for i = 1, . . . , n are i.i.d. and they are usually referred to as the complete or augmented data. Let y = ( y 1 , . . . , y n ). For a stipulated mixture complexity m ∈ N, let us denote by l c θ pm the log-likelihood of the complete data; i.e., On the other hand, the log-likelihood of the observed data x is given by It can be shown that the MLÊ θ p m = arg max θ pm ∈ pm l θ pm (x). (2.4) can be obtained by alternating between an expectation and maximization steps involving both the complete log-likelihood l c θ pm . This is precisely what the well-known EM-algorithm does. In the first step, the conditional expectation of l c θ pm ( y) given the observed data x is computed under the current parameter. Then, the obtained expression is maximized over the parameter space and the maximizer becomes the new parameter. These two steps are repeated until convergence. If s is the number of the iteration of the current E-step, then it is easy to show that this step is completed by computing the conditional expectation of the multinomial vectors z i given the observed data x. This yields for i = 1, . . . , n and j = 1, . . . , m Note that the maximizing mixing probabilities are easily obtained and are explicitly given in the s-th M-step by the expression . . , m, a numerical method might be required in case a a closed form is not possible. The optimization procedure then seeks to find at least the local maximum as finding the global maximum is not always possible. As noted in [50], the latter often occurs in the case of Gaussian mixtures with nonhomogeneous dispersions (unequal covariance matrices). Components that have either one observation, or several identical observations or several nearly-identical observations, result in the estimated covariance matrices that are singular, which causes the likelihood function to be unbounded. Gaussian mixtures with homogeneous components result in covariance matrices that are restricted in the parameter space and thus do not have this problem. For references on the EM-algorithm, see e.g. [23] and [48].
The description given above treats one given m, a candidate for the true mixture complexity. To obtain an estimator for m 0 , the true complexity, one can resort to maximizing a penalized version of the observed log-likelihood. This means that the log-likelihood will be augmented by a penalty term depending on the model complexity. Several widely used examples of this technique include Akaike Information Criterion (AIC) [2], Bayes Information Criterion (BIC) [62], Integrated Completed Likelihood (ICL) [10], Laplace-Empirical Criterion (LEC) [49], Normalized Entropy Criterion (NEC) [9] and many others [50]. These only differ in the form of the penalty function, and we will concentrate on the two criteria that have gained most popularity in practice: The Bayesian Information Criterion (BIC) and Integrated Completed Likelihood (ICL). While BIC is most widely used for performing model selection tasks, ICL is most frequently applied for solving clustering problems.
The general idea is to treat the task of choosing the number of components in the mixture as a model selection problem by considering a sequence of models where ν M m = p m is the number of independent parameters in the model andθ p m is the MLE of θ p m . The true complexity is then estimated by finding the integer m which maximizes BIC m . Given the discussion above, finding the number of components in the mixture that maximizes m → BIC p m is equivalent to choosing the mixture model with the greatest a posteriori probability. Some of the advantages of the BIC approach are that it is easy to implement, can be used for comparing non-nested models and was shown to be consistent for choosing the correct number of components in [40].
The ICL approach uses the log-likelihood of the complete data and replaces the unobserved labels z i j , 1 ≤ i ≤ n, 1 ≤ j ≤ m by their maximum a posteriori (MAP) estimator, that isẑ * i j = 1, ifẑ i j = arg max 1≤k≤mẑ ik 0, otherwise.
Thus, for the mixture model with m components The very useful relationship between BIC m and ICL m can be shown: It has been shown in [31] that in some cases (e.g. for the mixtures of Gaussians) evaluating the likelihood at the a maximum a posteriori (MAP) estimator instead of the MLE helps the EM algorithm to avoid singularities or degeneracies.
Regularization and variable selection techniques have also found their application in this setting. For example, [55] proposed an estimation technique for Gaussian mixtures in the context of a clustering problem, where the likelihood function is augmented by an L 1 -norm penalty term −λ m j=1 p k=1 |μ jk |, where μ jk is the k-th coordinate of the j-th mean vector, and derived a modification of an EM algorithm fitted for the purpose. The L 1 penalty can shrink some of the fitted means toward 0, thus leading to the most parsimonious model.
the E-step at the s-th iteration will update the probabilities given the current parameter vector θ (s−1)ẑ The M-step provides the following solutions: Further examples (for the mixtures of geometric and Poisson distributions) can be found in Appendix A in the supplementary materials.
Concluding this section it is necessary to point out that a great amount research has been carried out in this area, and multiple software applications have been developed for working with mixture models, in particular with the Gaussian mixture models that are most frequently used in practice. Most of the software is suited for model-based classification and in particular offering the opportunity to find the ML estimates via the EM algorithm. We refer interested practitioners to R packages Mclust [63] and mixtools [7] or the Matlab package MIXMOD [11].

Methods Based on the Hankel Matrices
The method of moments is generally considered to be less efficient when compared to maximum likelihood. Nonetheless, as justly argued in [46], there are situations where the method of moments reveals a nice mathematical structure. This is the case for the problem of estimating the true complexity of some finite mixture. As we will see below, the number of support points of a discrete mixing distribution with a finite number of jumps can be elegantly linked to whether the determinant of a special matrix of moments is equal to zero. Such a matrix is known under the name of a Hankel matrix.
We devote this section entirely to the estimation approaches based on the determinants of Hankel matrices of moments of the mixing distribution. The original method, with which we will start, was proposed in [21]. Additionally to the original approach we will describe a couple of its possible extensions. [21] motivated the method with a number of appealing features: 1. it gives consistent estimators under some mild conditions; 2. it requires no a priori upper bound on the unknown order of the mixture; 3. it comes with low computational time as it does not involve estimation of the mixture parameters.
Another attracting property, not mentioned by the authors, is that the method bears a universal character and can be applied to continuous distributions as well as discrete distributions with no modifications, provided that the moment generating function of the distribution exists.
For the reader to be able to appreciate the elegant argument standing behind the method, we shall recall next the key theoretical results furnishing its basis.

The Main Theoretical Results and Basic Approach
Recall that we have confined the present study to a one-dimensional case, where ⊆ R. For a given integer m ≥ 1 define the set In other words, the component c j is equal to the j-th moment of some distribution function G. For convenience, we will write c 2m = (c 1 , . . . , c 2m ) T for any given real numbers c j , j = 1, . . . , 2m. In [21] this set is defined more generally with nonnegative measure G.
For c 2m ∈ R 2m , the Hankel matrix associated with this vector is the (m+1)×(m+1) real symmetric matrix , denoted H (c 2m ) and given by Next we state the key result which links the true complexity of a finite mixture to the Hankel matrix of moments. See also Proposition 1 in [21]. Now we explain how the result above can be applied in the context of estimating the complexity of a finite mixture. Consider f 0 , a finite mixture of densities which belong to some family F and let G 0 be the associated discrete distribution function with true complexity m 0 . Then, Theorem 3.1 says that In other words, the correct order of the mixture is the first integer which sets the determinant to zero. But the theorem implies also that D m = 0 for all m ≥ m 0 . This characterizing feature of the true complexity is exploited to construct a sensible estimator. Indeed, assuming that it is possible based on the random sample (X 1 , . . . , X n ) to obtain a strongly consistent estimator of any j-th moment of G 0 ,ĉ j say, then the Hankel estimator of m 0 proposed in [21] is given bŷ {a m } m≥1 is a positive and strictly increasing sequence, and {l n } n≥1 a positive sequence satisfying lim n→∞ l n = 0 (we have omitted writing the subscript n in the notation of the estimators of the moments and determinants). Clearly, the term a m l n is acting as a penalty. Adding a penalty term to | D m | is necessary because otherwise minimizing of m → | D m | alone might yield an inconsistent estimator. In fact, strong consistency ofĉ j implies that | D m | is a strongly consistent estimator of the true value |D m | = D m (see our remark below). Since the latter is equal to 0 for all m ≥ m 0 , | D m | will be close to 0 for all m ≥ m 0 , which might result in choosing a value which is strictly larger than m 0 . Under some additional assumptions, consistency ofm n as defined above in (3.2) can be established as shown in Theorem 1 of [21]. We recall this result below.

Remark 3.1 Recall that
. Thus, D m seen as a multivariate real function of c 1 , . . . , c 2m (the components of c 2m ), is infinitely differentiable. Thus, ifĉ j is a strongly consistent estimator of c j for any integer j ≥ 1, then D m = det H (ĉ 2m ) is also a strongly consistent estimator of D m . Furthermore, a multivariate weak convergence ofĉ 2m toward c 2m as in the case where a multivariate Central Limit Theorem applies, the estimator D m will converge weakly to D m at a rate that is as fast as that ofĉ 2m . Typically, the estimatorsĉ j will result from considering some empirical estimators which we know to be asymptotically normal. Below, we will touch upon this point in some more detail.

Remark 3.2 In the light of Remark 3.1, the condition
Without going into the full proof of Theorem 3.2, let us give some intuition for the From the characterization if m 0 in (3.1), it follows that D m = 0 for m < m 0 implying that is assumed to be strictly increasing. Thus, we expect that as n → ∞ the minimum of the penalized criterion to be achieved at m 0 . The statement about consistency ofm n can be made more refined under additional regularity conditions. More precisely, suppose that for any integer m ≥ 1 there exist integrable functions ψ j and f j for j = 1, . . . , 2m such that the j-th moment of G 0 is given by T . Define the estimatorm n the same way as above witĥ We have the following theorem. See also Theorem 2 in [21].
Furthermore, assume that n 1/2 l n → ∞. Then, there exists a constant d > 0 and integer n 0 > 0 such that for all n ≥ n 0 The main argument in the proof uses judicious upper bounds on the probabilities P(m n ≤ m 0 ) and P(m n > m 0 ) based on concentration inequalities that involve the Cramer transform of the logarithm of the generating function of the centered random variables ψ j − E[ψ j (X )] for j ∈ {1, . . . , 2m} and m ≤ m 0 . Before commenting on the result itself, we would like to give some examples, which are relevant for the simulations section coming ahead. Example 2: Mixture of Gaussian distributions. Consider a finite mixture of Gaussian distributions with density If X ∼ f 0 , then X has the same distribution as Z + Y where Z ∼ N (0, 1) and Y ∼ G 0 with G 0 the mixing distribution with support points θ 1 , . . . , θ m 0 and mixing probabilities π 1 , . . . , π m 0 such that Y and Z are independent. Thus, for any j ≥ 1 where μ 0 = 1 and for an integer r ≥ 1 Thus, the vector of moments c 2m satisfies the triangular linear system c 2m In this case, we have Thus, for location mixtures of Gaussian distributions we have shown that More examples are available in Appendix A in the supplementary materials. Now we turn to commenting on Theorem 3.3. Although the result of that theorem seems to give an actual guarantee on the consistency ofm n , the exponential bound on the probability of being wrong about m 0 depends on n 0 and a constant d which are unknown. In case d is small and n 0 quite big, then consistency will not be observed for moderate and even big sample sizes: one would need an unrealistically huge number of observations to find the true complexity. Another problem is the estimation of the moments c j , j = 1, . . . , 2m for large values of m. Although the method does not require to put an upper bound on m while finding the minimum of | D m | + a m l n one has to choose some maximum admissible value for the mixture complexity. For large values of m the j-th moment c j can become very large. When this is combined with a low quality estimatorĉ j ,D m may be far away from 0, which is known to be the theoretical value for m ≥ m 0 . Such a phenomenon is illustrated using mixture of Gaussian distributions In Table 1 we give the first 8 theoretical moments c j of the mixing distribution and the mean value of their estimatesĉ j based on 100 replications for each of the sample sizes shown in the table. Table 2 gives the corresponding mean value ofD m as well as its penalized versions with a m = m and l n = log n/ √ n or l n = √ log n/ √ n for m ∈ {1, 2, 3, 4} computed on the basis of the same replications. It is clear from the values of Table 2 thatm n = 1 even for this very well-separated mixture and for the large sample size n = 10 4 .
The parametrization we chose implies that f 0 is a mixture of geometric distributions with success probability 0.7 and 0.2 respectively. In Table 3 one can see that the moments are accurately estimated for all sample sizes. However, the results of Table 4 indicate that the estimatorm n fails often to pick the correct mixture complexity, here m 0 = 2 for the penalties a m l n = m log n/ √ n and a m l n = m √ log n/ √ n. Our decision to take the penalty a m l n = m log n/ √ n is based on the recommendation made in [21]. The second penalty a m l n = m √ log n/ √ n was added in these simulations in order to compare the results obtained with the basic approach of [21] with the first modification we propose below and which is based on scaling the estimates of the determinants.  Table 4 The The inconsistency noted in these examples, despite the nice theoretical guarantees of convergence ofm n , are due to different reasons. In the Gaussian mixture, the penalty plays almost no role as | D m | dominates with values that are blowing up as we let m take larger values. For this reason, the estimator picks m = 1 in all cases. In the geometric mixture, the picture is completely reversed since the moments c j ∈ (0, 1) and hence get smaller for larger orders j. This causes | D m | to decrease with m. In this case, the penalty dominates and again m = 1 is often found as the minimizer of the penalized criterion.
The basic approach of [21] can suffer from serious underestimation (or overestimation) of the true mixture complexity even when the sample size is very large. Moreover, the question of how to choose the penalty term a m l n is not really settled in [21]. In fact, a penalty which would work for a certain family of distributions could perform miserably for another. A traditional approach here would be to resort to cross-validation to decide on an optimal choice for the penalty. Although this is a very important aspect of the problem, we choose not to pursue it here.

Modification of the Basic Approach Using Scaling
The discussion above reveals that while the estimator defined by [21] is quite appealing, it is very difficult to achieve consistency in practice. The main problem resides in that the method does not exploit any knowledge about the variability of | D m |. Without integrating the information about how these random variables behave stochastically (for n large enough), it would be almost impossible to say for example whether a small value obtained for | D m | can be seen as an indication that the true determinant is really equal 0. One way of circumventing the above issue is to use a rescaled version of this estimator. The starting point here is to use the already noted fact that the true determinant of the Hankel matrix of moments c 2m → D m is an infinitely differentiable function on R 2m . Thus, assuming that we can use the Central Limit Theorem to the vector of estimatorsĉ 2m , then for any fixed m > m 0 we get by applying the δ-method that Although the covariance matrix is unknown it can be estimated using resampling techniques. Here, we focus only on estimating the diagonal elements of , σ 2 1 , . . . , σ 2 m . By sampling B times with replacement from the original sample (X 1 , . . . , X n ) we obtain a new sample (X * 1 , . . . , X * n ) which can be used to compute the bootstrap determinants { D * 1 , . . . , D * m }. Repeating this procedure B times allows us to estimate σ j / √ n by computing the standard deviationσ * j of the bootstrap sample As the setting here is very standard, it follows that as n, B → ∞σ * j ≈ σ j , j = 1, . . . , m in probability. As mentioned in the previous section, the true order of the mixture can be assumed to be smaller than some given value m = m max ; i.e., the search of the minimizer of the penalized criterion will be performed in the set {1, . . . , m max }. Assuming that m max > m 0 , define the rescaled vector Thus, we redefine the estimatorm n aŝ We will not give a formal proof of consistency ofm n . The latter, however, can be intuitively seen to hold since it follows from the weak convergence in (3.6) that The proportions are computed on the basis of B = 500 and 100 independent replications Table 6 Proportion of the time the scaled Hankel estimator is equal to m 0 = 2 in the example of the finite mixture of geometric probability mass functions given in (3.5) n 100 1000 10000 The proportion is computed on the basis of B = 500 and 100 independent replications Gaussian vector with a covariance matrix having all its diagonal terms equal to 1. One the one hand, this implies that for any integer m ∈ {1, . . . , m 0 − 1} the probability that m is the location of the minimum should decrease as n → ∞. On the other hand, for m ≥ m 0 the penalty a m √ nl n becomes the dominating term. Since a m increases with m, the minimum is achieved at m 0 with increasing probability.
In the following, the examples considered above will be revisited using this modified approach to see to what extent the estimation accuracy is ameliorated. More specifically, we use the scaling approach described above to compute the proportion of times the alternative estimatorm n defined (3.7) is equal to the true complexity in 100 independent replications. In both examples, we have taken m max = 8. The number of bootstrap replications was taken to be B = 500. From Table 5 and 6 , we see how the results drastically improve with the scaling method for the sample sizes n = 1000, 10000 with 100% or close for the recovery of the true complexity. The improvement seems to be more pronounced with the choice of penalty m √ log n/ √ n. Thus, one conclusion that can be drawn here is that the method would greatly benefit from comparing the performance of different penalties. As mentioned above, such a comparison can be done using some cross-validation approach.

Modification of the Basic Approach Using Bootstrap
A specific feature of the Hankel matrix of moments methods discussed previously is the possibility to estimate the order of the mixture without estimating the parameters. However, it seems that there might be a high price to pay for avoiding this part: some of the essential features of the mixture may not be captured by the determinant alone, Table 7 Proportion of the time the modified Hankel estimator is equal to m 0 = 3 in the example of the finite mixture of Gaussian densities given in (3.4) n 100 1000 10000 Proportion of timesm n = 3 0 0 0 The proportions are computed on the basis of B = 500 and 100 independent replications which can lead to the wrong answer with a "bad" penalty, even for very large sample sizes. Furthermore, in many applications it may be desirable to obtain the estimator of the order of the mixture as well as the estimators of all the parameters. We describe here another modification that is suited for this purpose. It is in essence a sequential testing procedure in which some statistic computed from the data is compared with a critical value obtained e.g. by re-sampling from the assumed theoretical model. The said statistic can be taken to be either the determinant of the Hankel matrix D m as in the basic approach proposed by [21] or its rescaled version as described in the previous section. The idea is to replace minimizing the objective function in (3.2) or (3.7) by taking a reject/accept decision regarding whether the current value of m is equal to the true complexity. We describe this procedure only when the statistic is taken to be equal to D m since the modifications are rather obvious for the scaled version thereof: repeat the previous steps with m + 1.
The procedure described above is not new in the context of estimating a mixture complexity. In fact, a similar approach will be encountered below with the only difference that it is based either on some minimum distance estimation or likelihood ratio statistic (see Sects. 4 and 5 ). In a nutshell, one sequentially tests and declares as an estimate for m 0 the first value of m for which H m 0 is not rejected. In Table 7 and 8 we report the proportion of the time the sequential procedure described above gives the correct mixture complexity for the same finite Gaussian and geometric mixtures given in (3.4) and (3.5) respectively.
From the simulations results obtained in Table 7 and 8 we can see that this other modification of the original Hankel matrix method is less successful for the Gaussian mixture but still works well for the geometric one. This might be explained again by the large values of the higher-degree moments of the Gaussian distribution which  impact heavily the quality of estimating the determinants. This is not at all an issue with the geometric distribution whose moments are much easier to estimate.

Extension of the Hankel Matrix Approach using Neural Networks
The conclusions achieved on the basis of the simulation study, summarized in Sect. 6, stipulate that there is a need of search for a more reliable and universal mixture order estimation technique that would yield more precise estimates irrespectively of the underlying scenario. Obviously, the approaches we have already examined yield estimators which depend on the features involved in a non-linear fashion. To this extent we decided to turn our attention the popular statistical tool designed specifically for modelling nonlinearities between the sets of input and output variables -Neural Networks (NNs), in the hope that they might identify patterns and relationships that the other approaches cannot capture. For the past decade the amount of research carried out in the field of NNs has experienced exponential growth, and a great multitude of NN types and classes have been designed to successfully solve a wide range of problems. The simplest of these tasks like image labelling or pattern recognition are usually solved by feed-forward network architectures such as the multi-layer perceptron (MLP), convolutional neural network (CNN) or radial basis function network (RBFN). More sophisticated tasks such as speech recognition or text translation require more complex interactions between the layers of the network, which are implemented in such architectures as long-short-term memory (LSTM).
At this stage of our research we are not aiming at estimating the whole mixture model (finding the optimal complexity as well as all its parameters) but only pursue the goal of identifying the number of subgroups in the population based on the observed data. Thus, when cast into the NN framework, the problem of estimating the number of components in a mixture can be viewed as a supervised multiclass classification problem and the relevant questions that need to be addressed are • discovering the most informative features to be fed into the NN and • devising an adequate design for the training set; • proposing the optimal NN architecture; • choosing the learning algorithm; • determine whether a universal architecture for multiple families of distributions can be found; • understanding whether using NNs is beneficial compared to the other methods. It seems natural to start the search for a well-suited network by browsing first thorough the simplest class of the NNs: the multilayer perceptron (MLP). Despite of their relative simplicity, networks with just two layers can approximate any continuous functional mapping [12]. One of the simplest possible architectures of the considered model with only two hidden layers is depicted in Fig. 1.
The input features x 1 , . . . , x d , d ≥ 1, the first and the second hidden layers consist of u and r neurons respectively, the weights vector ω is learned by optimizing the loss function J (ω = (ω (1) 1,1 , . . . , ω (3) 6,r ) T ), which is taken to be the cross entropy function, which is most often used in multiclass classification tasks such as ours: where z k = 1, if m 0 = k, with k = 1, 2, 3, 4, 5, 6, 0, otherwise are known labels for each of the generated vector of features in the training sample.
To this end, the choice over the optimal configurations is restricted to deciding on the number of hidden layers in the network, the number of neurons in each layer and the corresponding activation functions, the loss function and the learning algorithm.
For estimating the order of a mixture using a NN, one needs an appropriate assumption on the possible maximal number of components. We take the maximum number of components to be 6 for our task.
Recall from Sect. 3 that Hankel's criterion is backed by elegant, orderly statistical theory, however the method's performance turns out to be rather poor in practice. We use Hankel matrix determinants as inputs to the MLP to try to improve the estimation results by exploiting the information concentrated in the determinants without resorting to the use of any penalty function. Our experience reveals that using sequences of the first 6 Hankel matrix determinants of the mixing distribution as inputs leads to improved results when compared to other tested options. For this reason we regard this approach as an extension to the original Hankel technique.
Using the sequence of the Hankel determinants as inputs produced resulted in high performance for the Geometric mixtures, but showed poorer performance for the Poisson mixture due to the fact that the determinants for the Poisson mixtures tend to blow up while those for the Geometric mixtures stay bounded within a [0, 1] interval. For that reason the inputs for the Poisson mixture had to be modified. One of the modifications led to good performance was the relative changes in the absolute values of the Hankel determinants: To ensure the variety of the training examples, the characteristics and structure of the data that is used for prediction should be scrutinized and taken into account. The training set should be designed to be as representative as possible of the data of interest. The following procedure can serve as a useful example. For the considered mixtures distributions the parameters for each mixture component is chosen without replacement from a grid with a pre-specified step to insure that these are distinct. A grid on [0.05, 1] with step 0.05 for the geometric distribution should do well in most of the applications. A grid on [1,20] with step 1 for the Poisson can be taken if the expected rate of occurrences is believed not to exceed 20 by much and the parameters of subpopulations are separated by at least 1 unit. The step value can be reduced if more precision is desired. The mixing proportions can be taken on a gird with an appropriate step size in a similar way, in this case replacement is allowed and the generated results should be normalized. The mixtures with the parameters obtained in this way are then used for generating samples. It seems to be useful to enrich the training data set by simulating several times from each of the resulting mixtures to account for possible variation in the sample populations during the training.
The 6 output neurons of the output layer is further processed by softmax activation function σ : R 6 −→ [0, 1] 6 , which ensures that the estimated class probabilities live between 0 and 1 and sum up to 1: whereq k us the resulting value of the k-th output neuron. Whenever optimizing the loss function, the value of the learning rate becomes of importance: when too small, the weights of the NN are hardly updated, and much time is needed for the network to find a solution; when too large, the weights are updated in large increments, and overshooting the optimum becomes highly probable. In the case of mixture order estimation the value of the learning rate is influenced by the values of the parameters of the mixtures in the training set: an efficient learning rate for the geometric mixtures, where the parameters lie within (0, 1] and thus higher precision is required to separate the components, will be smaller than for Poisson or Gaussian mixtures where the parameters can range from 0 to 20. The output of the MLP is a vector of estimates of the class probabilities, that is, the probabilities of an observation (represented by either a sample from a mixture distribution or a vector of alternative relevant features) belonging to each of possible 6 classes:p k , k ∈ {1, 2, 3, 4, 5, 6}. The predicted number of components in the mixture is taken to be the class with the highest estimated probability: The search for a successful model (done using KerasTuner library [54]) requires examination of a large parameter space even for a simple network such as a MLP. Combinations of several hyperparameters of the NN were kept track of in order to identify the optimal NN architecture: • activation function: relu, tanh, sigmoid • number of layers: 1, 2, . . . , 9,10 • neurons in each hidden layer: 10, 25, 40, . . . , 280, 295, 310 • dropout layer after the last hidden layer: dropout rate between 0 and 0.1 • learning rate for the optimization algorithm: 10 −2 , 10 −3 , 10 −4 A set of 10000 samples was used for training the NN, the motivation being that taking the mo-ments and determinants as the features requires quality estimation thereof. Therefore, while a sample of this size is rarely available in real-world datasets, the emphasis was placed on finding a neural network that would perform well if the estimates are good. In practice, a neural net-work trained with 10000 samples still predicts well when the test sample size is much smaller. Unfortunately, we were not able to find a single MLP specification that would work equally well for all considered families of distributions -Gaussian, geometric and Poisson. Table 9 presents three different MLP configurations for Gaussian, geometric and Poisson mixtures that achieved satisfactory performance in our simulations. The predicted class probabilities as well as prediction accuracy on a number of test cases for networks with the denoted specifications can be found in Sect. 6.

General Setting
The estimation techniques discussed in this section are mainly based on the works [69,74,75]. Additional relevant references will be mentioned below. In a nutshell, these techniques use the minimal distance between a consistent estimator of f 0 ,f n say, and an estimator of the latter in the class of finite mixtures with m components m ≥ 1. Using the notation of Sect. 2, this means that the estimator of the true m 0 will be based on the projection off n on the class F m , m ≥ 1 in a sense that will be determined. As noted above, these classes are nested; i.e., F m ⊂ F m+1 . Thus, the basic approach stops at the first m where the projection on F m+1 does not bring a substantial improvement over the projection on F m . One main feature of the methods investigated in this section is that one performs estimation of the parameters of the mixture as well as the mixture complexity.
In the above mentioned papers, the projection of F m makes use of the Hellinger or the L 2 distances. The estimation method suggested in [74] is built upon a model selection procedure by sequentially fitting the nested mixture models. Thus, the method is reminiscent of the sequential testing approach already encountered above for the second modified version of the determinant of the Hankel matrix of moments. The procedure allows at each iteration to search over a higher class by adding one component to the mixture and find the best model within each class until adding another component brings no more benefit in the sense that it decreases the objective loss function by an amount smaller than a specified tolerance level. In the sequel, we consider only the situation where the dominating measure, μ, is either the counting or Lebesgue measure. In the first case, the nonparametric estimator of f 0 is the empirical probability mass function given bŷ In the second one where X is absolutely continuous, we consider a kernel density estimator with fixed or random bandwidth c n : where K is some standard kernel function.

The Minimum Distance Estimator: The Basic Approach
In the following, let F denote the class of densities that are mixed. Also, let D be either the Hellinger or L 2 distance, that is for two densities f and g with respect to μ we have either Recall the following notation from Sect. 2: For a given For a given density f , we now consider the functional provided that the minimum exists. Here, θ D p m ( f ) denotes the set of all minimizers as uniqueness of the solution is not guaranteed. Note that our notation is different from the one used in [8], [74], [75] and [69]. In case f =f n , we have the following definition: for a given m and the non-parametric estimatorf n defined above in (4.1) or (4.2), the minimum distance estimator with respect to D of θ p m is defined aŝ provided that a minimizer exists. Proving existence of a minimizer θ D p m ( f ) for some given density f requires careful argumentation under some regularity conditions. When D is the Hellinger distance, Theorem 1 in [8] gives a proof of this existence under the condition that t p m → f t pm (x) is continuous for almost every x ∈ X , that the mixture is identifiable and the parameter space p m = S m−1 × m is compact. Also, the same theorem proves uniqueness of θ D p m ( f ) in case f is itself a finite mixture. In order words, if f = f θ pm , then θ D p m ( f ) = θ p m . This can be easily seen as an immediate consequence of identifiability. When D is the L 2 distance and μ is the counting measure on the set of non-negative integers, then [69] show a similar theorem while relaxing the condition of compactness on the parameter space. The main building block in the proof is to show that the mapping is continuous. In the discrete setting considered in [69], a proof of the continuity property can be based on a slightly different argument. Indeed, if t (k) p m is a sequence converging to t p m as k → ∞, then by Minkowski's inequality (also used in page 4252 of [69]) using that a pmf is always bounded by 1.
The latter sum converges to 0 by continuity of t p m → f t pm (x) and application of the Sheffé's Theorem. For existence of a minimizer when D is the L 2 -distance, [69] makes the assumption that for the pmf f to be projected, for any m there exist some compact C (which depends on m but we omit writing this dependence explicitly), and θ * p m such that inf g∈ pm \C Such an assumption is not needed in case p m is itself compact. Also, the compact C is rather abstract and one only needs to exhibit its existence in some way. It is clear that even when θ D p m ( f ) is not a singleton, we have f θ pm = f θ pm a.e. for two minimizers θ p m , θ p m ∈ θ D p m ( f ). When f =f n , we will denote the corresponding density f θ pm (or f θ pm ) by f D m . Note that we have omitted the subscript n, and replaced p m by m for the sake of a lighter notation. By definition of θ D p m ( f ) we have The roles of f D m and f D m will become clear below. Although our notation for those projections is different from the one used in the aforementioned papers, our choice is driven by our desire to maintain some notational coherence throughout this survey. Now, we describe how the basic approach works with the minimum distance estimators. The estimation procedure as outlined in [74] is much inspired by the work of [37]. The latter paper is however mainly focused on estimating the true complexity of a finite mixture of Gaussian distributions using the Kullback-Leibler divergence instead of the Hellinger (or L 2 ) distance. There are two starting points of the basic approach. The first one is to note that the true mixture complexity m 0 satisfies While the identity in (4.6) is a direct consequence of identifiability, the one in (4.7) is less obvious. Note that this identity is proved if we show that for m ≤ m 0 − 1, . A proof of this fact when D is the Hellinger distance can be found in p. 1485 of [74], where an auxiliary lemma (Lemma A.4) was used. For the case where D is the L 2 distance, and for the sake of completeness, we give a proof in Appendix C in the supplementary materials.
Based on the discussion above, it seems natural to search for the first m which minimizes some empirical version of the distance . The second one is to recall once more the inclusion F m ⊂ F m+1 . This implies that The inequality above means that without penalization it is in principle impossible to find a finite order which can be taken as an estimator of m 0 . This overfitting is accounted for by adding a penalty term which is proportional to the number of parameters in the mixture model. This yields the following criterion for some chosen sequences {v m } m and {b n } n such that the former is increasing and the latter satisfies lim n→∞ b n = 0. Note that in the works [74], [75] and [69], b n v m /n is taken instead. Now, mimicking the property in (4.7) gives rise to the following definition of the minimum distance estimator If this minimum does not exist, thenm n = ∞. The term α n,m can be seen as a threshold so that an integer m is declared to be the estimator when the projection of f n on the class F m+1 yields an insignificant change in comparison with its projection on the previous class F m .
Strong consistency of the minimum distance estimator defined above result was stated in Theorem 1 in [75] for the Hellinger distance and in the Consistency Theorem Section in [69]. In the former, the nonparametric estimatorf n is taken to be a kernel estimator, as defined above in (4.2), with bandwidth c n such that lim n→∞ c n + (nc n ) −1 = 0. In the latter,f n is the empirical probability mass function as given in (4.1). Now, the only condition assumed on α n,m in these convergence theorems is that α n,m → 0 as n → ∞. This is of course guaranteed by the fact that lim n→∞ b n = 0. However, we believe that this statement is not accurate as is, since the penalty needs to depend on the rate of convergence off n to the true density f 0 . This rate of convergence is known to depend on the smoothness of f 0 and the kernel K . In Appendix C in the supplementary materials, we explain why the condition lim n→∞ α n,m = 0 is not enough.

Modification of the Basic Approach Via Bootstrap
In this section we propose a modification of the basic approach based on minimal distance between the projection of a non-parametric estimator on the class of mcomponent mixtures and this estimator augmented with some given threshold. We resort in this modification to a parametric bootstrap procedure in order to avoid the bad choice of a threshold. For a given integer m ≥ 1, consider the hypothesis testing as in (3.8) Ifq B,α/2 andq B,1−α/2 denote the empirical α/2 and (1 − α/2)-quantiles based on the bootstrap sample and the current candidate m is replaced by m + 1. We take as our estimator for m 0 the first m for which the null hypothesis is not rejected. Here, we report simulations results for the two examples for finite mixtures considered above; see (3.4) and (3.5). In these simulations, B = 500 and the sample sizes considered are n = 100, 1000, 10000. The performance of the method proposed in this section is assessed through the proportion of times the estimator is equal to the true complexity (3 and 2 respectively). We used Hellinger distance for the Gaussian mixture and L 2 distance for the geometric mixture.

Sequential Likelihood Ratio Tests with Bootstrap
Likelihood-based methods play a central role in statistical inference. Among these methods the likelihood ratio test (LRT) is one of the most widely used in practice. The LRT has a simple interpretation and enjoys the property of being invariant under re-parametrization; see for example [41]. Furthermore, whenever certain regularity conditions are satisfied Wilks' theorem [71] says that, under the null hypothesis, it converges weakly to a chi-square distribution as the sample size grows to ∞. The degrees of freedom of the limiting chi-square distribution are determined by the difference of the dimension of the whole parameter space and that of the null space. However, and as already mentioned above, the regularity conditions required for the asymptotic theory to hold are not satisfied in the mixture problem. The issue is that it is always possible to write a mixture model with m components as a model with m + 1 components (or more) by setting some of the mixing probabilities to 0. Hence, the true parameter under the null hypothesis lies on the boundary of the alternative space. To better explain the problem, let us consider the example of testing whether a random sample was generated from a (single) Gaussian distribution N (θ, 1) versus a 2-component Gaussian mixture, where each of the components has variance equal to 1. If f 0 denotes the true density, then the goal is to test for some θ ∈ R and θ 1 = θ 2 ∈ R and π 1 ∈ (0, 1). If we define here the likelihood ratio as the ratio of the maximum likelihoods over the null and the whole space, then the maximization needs to be done on R (the null space) and (equal to the whole parameter space; i.e., the union of the null and alternative spaces). Thus, a density under H 0 lies on the boundary in the sense that f 0 results from either setting the mixing probability π 1 to 1 or 0 and taking θ 1 = θ (or θ 2 = θ ), or letting θ 1 = θ 2 while π 1 is arbitrary. In the classical setting, one of the main arguments which leads to the chi-square distribution as the weak limit is the use of Taylor expansion up to the second order of the log-likelihood at the global MLE around the true parameter. We can refer here to the proof of Theorem 22 in [26] for well-explained and rigorous arguments. Things then work because the true parameter is an interior point and has a unique representation as an element in the whole parameter space. In the mixture model setting, this argument does not work because, as it can be seen from the example above, the true parameter has different representations under the null hypothesis where it is on the boundary. This non-standard situation has triggered a strong interest for either computational or theoretical investigation of the limit distribution of the LR statistic under the null hypothesis. In this context, we can refer to [1], [64], [13], and [47]. For a nice review of the papers on this subject, we can refer to Section 6.5 in [49] among others. The main message that one can take from these works is that when the limit distribution of the LRT can be simply described, it is a mixture of chi-square distributions. In more complicated cases, this limit distribution is given more abstractly by sup s∈S max(0, Y s ) 2 where Y is some well-defined centered Gaussian process and S is some suitable set. See for example [47] where it is shown that S is the ensemble of cluster points of some generalized score.
In this section we review the sequential testing procedure based on LRT and use of resampling for approximating the distribution of the test statistic under the null hypothesis. This method was proposed in [38] for estimating the true complexity of a finite mixture of Poisson. Although the focus there was put on that family, the approach can certainly be extended to other distributions. Consider again the hypothesis testing problem where L X denotes the likelihood function based on the sample X = (X 1 , . . . , X n ) from the unknown mixture, and p m = m(d + 1) − 1. As in [38], consider the loglikelihood ratio statistic A mixture with smaller number of components will be rejected in favor of the larger model whenever the log-likelihood ratio statistic is large. Otherwise, the mixture will  In view of the issues related with deriving the asymptotic distribution of the loglikelihood ratio statistic, one can resort to a parametric bootstrap. As this is already done above, the description of the procedure is omitted.
We give the results of this procedure for the two finite mixtures with m 0 = 3 and 2 already considered above.

Simulation Results
In this section we describe the simulation results obtained for sample sizes n ∈ {50, 100, 500, 1000, 5000, 10000} using the procedures described in the previous sections for finite mixtures of Gaussian, geometric and Poisson distributions with m 0 ∈ {2, 3, 4}. Throughout this study, unless explicitly stated otherwise, the standard deviations for Gaussian mixture components are assumed to be 1. For the minimum distance-based methods, we follow the recommendations made in the relevant papers [75] and [69] and use the following two penalty functions α n,m , with m denoting the stipulated mixture order and n the sample size: • the penalty based on the Akaike Information Criterion (AIC) α n,m = 0.6 n log m+1 m for the L 2 distance and α n,m = 2 n for the Hellinger distance • the penalty based on the Schwarz Bayesian Criterion (SBC) α n,m = 0.6 log n n log m+1 m for the L 2 distance and α n,m = log n n for the Hellinger distance.
For all simulations, we used 500 replications to compute the frequencies of an exact recovery of the true complexity. These are displayed in Tables 28, 29, 30, 31, 32, 33 and 34 to be found in Appendix B in the supplementary materials.
Clearly, the performance of the NN's cannot be directly compared to the performance of the other methods we have considered. One of the reasons being that the NN framework is very much different from the other approaches in terms of training and testing procedures applied. The neural network requires to encounter as many different mixtures as possible to be able to learn the latent structure, while the other techniques (except for the Hankel approach) are seeking the best fit for a given sample from the target mixture. We still report the NN results along with the performances of the other methods with a few reservations so that the reader could put together the full picture. Several points need to be borne in mind when reading and interpreting the NN performance: • the values reported are the predicted class probabilities that reflect how sure the network is that an observation (one of the tested mixtures) belongs to each of the represented classes; • bin ≥ 5 can be slightly misleading in this case as it reflects the summed probability for classes with 5 components and 6 components and in a few instances the probability corresponding to one of the classes is larger then two of the probabilities for 5 and 6 when considered separately but less, when they are so combined; still the results as given provide an insight into how well the NN can tell the classes apart; • sample size of 10000 was used for computing the input features for the NN thus the results are reported in the respective table cell; • 10000 samples were used to train the NN; • accuracy measure that shows in how many instances the NN was able to correctly predict the true number of components cases will be reported separately.
The accuracy measure was computed by generating 100 different samples from one of the selected mixtures (these did not occur in the training set), and the percentage of times the network gave out the correct prediction was computed.
It follows from the simulations that the bootstrap modification of the Hankel matrix determinants method shows reasonable results for the geometric and Poisson mixtures but fails to produce accurate estimates for all 2-component and 3-component Gaussian mixtures, persistently underestimating the number of components when compared to the truth. This underestimation is likely to be the result of the "exploding"moments and determinants issue inherent in the Hankel matrices of the mixing distribution for the finite normal mixtures. The issue was covered in Sect. 3, and the simulation study results illustrate it once again.
For a well separated 2-component Gaussian mixture all methods (except for the Hankel matrix determinants approach with bootstrap) perform well even for sample size as small as n = 50. The BIC and ICL methods demonstrate high performance in this setting (estimating the number of components correctly in more than 95% of the cases for small sample sizes of n = 50, 100 and 100% for larger ones) and so do the minimum Hellinger distance methods. For the estimation of mixture complexity via the minimum Hellinger distance procedures (the L 2 distance is not applicable in the case of continuous distributions) we used a KDE with bandwidth of 0.5. The Hankel matrix determinants modification with bootstrap and scaling show slightly inferior results for small and average sample sizes when compared to BIC, ICL and the minimal Hellinger distance methods but performs well for sample sizes of n ≥ 5000. The LRT approach for the available sample sizes of n ∈ {50, 100, 500, 1000} shows more modest results in this setting although the relative frequency of the correctly estimated cases is still high. For a well-separated Gaussian mixture with a low mixing proportion for one of the components the LRT method works well irrespective of the sample size estimating the number of components correctly in more than 90% of the instances for n ≥ 100 and outperforming many other methods (except for BIC and ICL) for small sample sizes. The performance of all techniques drops significantly for not well-separated mixtures (e. g. a 2-component Gaussian mixture with overlapping regions). The NN extension of the Hankel method demonstrates good results in all 3 cases for the 2component Gaussian mixtures even when the modes are located close to each other, which is manifested in the estimated class probabilities as well as in high accuracy rates of 100%, 96% and 100% respectively.
The BIC and ICL methods that proved to be supreme when compared to other methods in many scenarios for the Gaussian mixtures are not performing as well for the mixtures of geometric distributions. For the 2-component mixtures the simulation results indicate that the minimum Hellinger distance with bootstrap and LRT techniques tend to outperform all other tested methods. The penalized Hellinger and L 2 -distance methods perform well whenever the mixture is well-separated (the weights do not have to be well-balanced however) and n ≥ 500. For small and medium sample sizes the minimum distance-based methods with penalties show rather poor results. The techniques using the Hankel matrix determinants also enjoy high performance for well-separated mixtures when the mixing proportions are similar or not of observations is large, identifying the number of components correctly in 95% of the cases. The NNs also demonstrate less confidence when applied to the geometric mixtures as indicated by a higher level of spread in the estimated class probabilities although the accuracy rate of 59% can be considered a rather decent success measure for the challenging first scenario. The accuracy amounts to 100% for a less demanding case with the well-separated mixture.
Generally, for n ≥ 5000 and more observations all methods with an exception of the minimal L 2 distance approach using the SBC-based penalty allow for the correct estimation in at least 90% of the cases.
For well-separated 2-component Poisson mixtures all methods perform well even for very small sample sizes. The BIC, ICL and the two penalized MHD methods often demonstrate accuracy which is close to absolute. The minimum distance-based approaches with the AIC-based term seem to work better than the minimum distancebased approaches with the SBC-based penalty whenever the mixtures are either not well separated or when one of the mixtures is scarcely represented. Whenever the mixture is less challenging, the approaches involving the SBC-based penalty tends to lead to higher accuracy than those involving the AIC penalty. Also, choosing the Hellinger distance for the minimal distance techniques seems to achieve slightly better performance than using the L 2 distance (whenever the components are not too close). In the case of well-separated 2-component Poisson mixtures the scaled bootstrap version of the Hankel matrix approach seems to outperform the bootstrap modification. For a well-separated mixture with a low mixing proportion of one of the components the LRT and the scaled version of the Hankel matrix determinants techniques seem to be an appropriate choice. Approaches using the Minimum L 2 distance perform rather poorly for all sample sizes in this setting while all Hellinger distance-based approaches provide good estimates when the number of observations is large. The NNs seem to be able to learn well the 2-component Poisson mixtures. The accuracy rates for 3 out of 4 scenarios are 100%. The only exception being the mixture with very closely located components, where the NN fails.
For a well-separated 3-component Gaussian mixture with similar mixing proportions the BIC and ICL methods show outstanding performance even for small sample sizes. All methods (except for the Hankel matrix determinants with bootstrap) perform relatively poorly for small sample sizes and very well for sample size of n ≥ 5000 where their accuracy achieves 95%. The Hankel matrix method with bootstrap and scaling shows again slightly worse results for small and average sample sizes but performs well for sample sizes of n ≥ 5000. The minimum Hellinger distance approach with AIC-based term shows better performance than other estimators from the same group for average sample sizes. The LRT method delivers poor results for small sample sizes but achieves mor tha 95% accuracy already for n = 500. For a well-separated mixture with a low mixing proportion for one of the components the picture is similar, with a slightly lower performance. The BIC method outperforms the other methods, although giving correct predictions in more than 80% of the cases only for n ≥ 5000. The LRT method seems to outperform all other methods (except for the BIC) for small sample sizes, providing the correct estimate in almost 70% of the instances for n = 500. In general the results for almost all methods are rather poor for small sample sizes, improving as the number of observations grows. For average sample sizes the minimum Hellinger distance with AIC-based term shows a significant improvement, estimating the mixture complexity correctly in more than 70% of the cases, and all Hellinger distance-based methods perform well for sample sizes n ≥ 5000. The NN demonstrates again good results in all 3 cases for the 3-component Gaussian mixtures achieving accuracy rates of 90%, 96% and 100% respectively for the given cases. For a not well-separated mixture all methods perform quite poorly, being able to identify only 2 components out of 3 in most of the instances. The LRT method for the available sample sizes shows estimation results which are only marginally better than those achieved by the other methods. The BIC approach also performs poorly in this case, often either underestimating or overestimating the number of components in the mixture, but still shows better results when compared to other methods.
Geometric mixtures with 3 components seem to be a challenging task for all methods considered in the survey. The minimum Hellinger distance approach with bootstrap and the LRT technique are still able to show better results for small and average sample sizes when compared to other methods whenever the mixtures are well-separated. It is likely that these methods could show better performance for large sample sizes. For the first 3-components mixture, which is not well-separated, none of the methods estimate the complexity correctly, systematically underestimating it. In this scenario the methods penalized with AIC-based penalty once again show better performance than their counterparts using the SBC-based thresholds. The group of distance-based approaches as well as the cluster of Hankel matrix techniques also perform poorly in this setting even when the sample size is large and the mixture is well-separated. The 3-component mixtures seem to be challenging also for the NNs, causing the accuracy rate for the two scenarios to drop to 65% and 44% respectively.
In the case of well-separated 3-component Poisson mixtures the BIC method achieves absolute performance for n ≥ 5000, while the LRT technique outperforms all other tested methods for small sample sizes. One could suppose that it would be very effective for larger sample sizes. The ICL persistently underestimates the number of components in the 3-component mixture even when the components are far apart. Whenever the Poisson mixture is well-separated the modification of the Hankel matrix approach with bootstrap tends to be more accurate than the Hankel matrix approach with bootstrap and scaling. All methods relying on the minimum Hellinger distance estimate the complexity correctly in more than 90% of the instances for n ≥ 1000 and the LRT techniques perform well for average sample sizes. The bootstrap modification of the Hankel matrix determinants approach and the minimal L 2 distance-based methods also perform well whenever the mixture is well-separated and the mixing proportions are approximately the same for all components, with no component dominating over the others provided the number of observations is large enough. For average sample sizes and "unbalanced" mixtures the results of these methods are rather poor when compared to the previously listed approaches. The NNs achieves accuracy rates of 83% and 76% for the two well-separated mixture, but is not able to correctly predict the order of the mixture with very closely located components.
Poisson mixtures with 4 components seem to be quite a challenge for all methods, and a large number of observations is needed to obtain decent performance.

Applications to Real Data
In this Section we demonstrate how the approaches discussed in the previous sections perform on real data. Wherever applicable the size of the bootstrap sample size is taken to be 1000, for the LRT-based approach, the observed LRT statistic is always compared with the 95%-quantile, for other methods using bootstrap (the methods relying on the Hankel matrix determinants and the minimum distance-based procedures), 2.5%-and 97.5%-quantiles are used.
We will begin the analysis by considering the Old Faithful Geyser Data, first published in [4]. The data comprises waiting times between eruptions and the duration of the eruption for the Old Faithful geyser in Yellowstone National Park, Wyoming, USA. The waiting times between the eruptions of the Old Faithful geyser are assumed to be well described by a finite Gaussian mixture model.
Unfortunately, the techniques based on the Hankel matrix determinants do not appear to be an appropriate choice in this setting. To make use of the translation nonparametric Hankel method approach (as outlined in the Section 3.1 of [21]) one should make sure that the assumption of equal standard deviations throughout the mixture components holds. This assumption does not seem to hold for the Old Faithful data. The NN extension cannot be used for these data either as for the NN training we assumed a known variance of 1 for all components in the mixture, which is not the case here.
The BIC and ICL method applied to the Old Faithful data result in different estimates of the optimal number of components. The maximal BIC values corresponds to 3  Table 14.
Implementing the minimum Hellinger distance method with an automatically selected bandwidth for the KDE and AIC-based penalty one obtains that the optimal mixture should have 2 components. The resulting estimated 2-component mixture (in green) as well as the 2 components (in dark red) are plotted in Fig. 2 along with the empirical distribution of the waiting times.
Changing the penalty term to the SBC-based from the AIC-based yields a slightly different parameter vector estimate, however the choice of the number of components remains the same. The differences of the squared Hellinger distance values and the AIC/SBC-based thresholds can be found in Table 15. The minimum Hellinger distance method with bootstrap and the LRT approach also yield a 2-component mixture.
The minimum distance methods and the LRT attain slightly different parameter estimates (the estimated parameters for BIC and ICL approaches are similar to those of the LRT as in all these cases the MLE is used). The results have been gathered into a single table (Table 16) for the ease of comparison. We now consider the data taken from the 1952 Annual Report of a pension fund that contains the information on the number of children of 4075 widows entitled to the fund support. The dataset first appeared in [67]. The data do not appear to be simply a random sample from a Poisson distribution as the number of zeros (widows with no children) appears to be too large. This issue was treated in [67] by fitting a mixture of two processes, one of which is a Dirac distribution at 0 while the other follows a Poisson distribution. We attempted to fit a mixture of Poisson distributions to the data using several of the discussed approaches to verify the above mentioned population heterogeneity assumption.
The BIC and the ICL methods estimate 2 and 1 components respectively for this data set, the values are given in Table 17.
Non-parametric non-scaled and scaled Hankel matrix approaches with respective penalty terms m log(n) √ n and m log(n) yield the estimated number of components in the Poisson mixture equal to 2, as can be seen from Table 18.
The parametric non-scaled Hankel matrix determinants approach with bootstrap as well as the corresponding scaled version yield a 2-component mixture of Poisson distributions with the estimated by the maximum likelihood parameters (0.66, 0.34, 0.03, 1.11). This agrees with the population's heterogeneity hypothesis  found in [67]. The estimated Poisson mixture along with the empirical distribution are shown in Fig. 3. The estimated parameters for the minimum Hellinger and L 2 distance methods with AIC-based and SBC-based penalties and bootstrap as well as the LRT approach also yield the same number of components and very similar parameter estimates. Table 19 displays the differences of the squared Hellinger and L 2 distances along with the corresponding thresholds.
The NN predicts 2 components with very high probability, the estimated class probabilities output by the NN being as reported in Table 20.
Thus all methods (except for ICL) applied to these data agree on the optimal number of components in the mixture, also confirming the hypothesis, identifying the optimal number of components as 2. Bold are those obtained for the true complexity of the mixture The dataset used for fitting a mixture of geometric distributions is the Shakespeare dataset, analyzed in the seminal work [24], which comprises counts of the number of times certain words that William Shakespeare used in his writings. The data is designed as follows: the number of times Shakespeare used a word only once is 14376, the number of times the same word occurred exactly 10 times in his writings is 363 an so on. The goal set in [24] was to use the observed frequencies words, to estimate the unobserved number of words that Shakespeare knew but did not use in his writings. This problem is known under the name of "species richness"and can be solved using a variety of approaches. One of the approaches was considered in [6], where the theoretical rationale for using a mixture of geometric distributions in such a setting is laid out.
The BIC and ICL methods both select 2 as the optimal complexity for the Shakespeare dataset as can be deduced from Table 21 containing the BIC and ICL values for number of components 1, . . . , 5.
Both non-parametric non-scaled and scaled Hankel matrix techniques estimated the number of components as 2, as follows from Table 22.
The parametric Hankel matrix determinants approach with bootstrap applied to the Shakespeare data yields the estimated number of components in the mixture of geometric distributions equal to 2. The same number of components is obtained when the scaled version of the parametric Hankel matrix determinants approach with bootstrap is used.
The estimated mixture for the minimum Hellinger and L 2 distance methods with AIC-based and SBC-based penalties and bootstrap suggest 3 components unlike the Hankel matrix determinants procedures, yielding slightly different estimates (0.4148, 0.3758, 0.2094, 0.8406, 0.2902, 0.0490) (for Hellinger) and (0.3839, 0.3870, 0.2291, 0.8622, 0.3203, 0.0523) (for L 2 ). The squared differences for Hellinger and L 2 distances and the BIC/AIC-based thresholds can be found in Table 23.
The estimated mixture of geometric distributions using the Hellinger distance approach with bootstrap, all of its components and the empirical distribution are plotted in Fig. 4.
The LRT approach results in a 3-component mixture with the estimated parameter vector (0.4270, 0.3744, 0.1986, 0.8311, 0.2741, 0.0446).  The estimated parameters for the minimum Hellinger and L 2 distance methods with AIC-based and SBC-based penalties and bootstrap as well as the LRT approach also yield the same number of components and very similar parameter estimates.
The NN predicts the maximal possible number of components for the Shakespeare data, which is 6, with the class probabilities designated in Table 24: Thus for the Shakespeare data various methods do not agree on the optimal number of components, which ranges from 2 components to 6, also providing different estimates for the model parameters. The final model choice in this case is left to the researcher who might have an insight on which number of component makes the most sense (e.g. according to the parts of speech different words belong to). Table 25 summarizes the estimates for all reviewed methods and all datasets.

Other Work on Mixture Complexity Estimation
The current survey certainly does not cover all the existing methods for estimating the complexity of a finite mixture.
Before mentioning some other interesting references, we would like to draw the reader's attention to the seminal work of Bruce Lindsay which, among other things, brought a very novel way of viewing the nonparametric maximum likelihood estimator (NPMLE) of a general mixture distribution; see [44] and [45]. The novelty resides in considering this estimator from a geometric perspective. One of the most important results which derive from this is that, under some simple conditions, the NPMLE of a general mixture is a finite mixture with complexity not exceeding the number of distinct observations. Not surprisingly, all the papers on the methods reviewed and implemented in this survey refer to one of Lindsay's work on mixture models. This continues to hold true for the other approaches we would like now to mention and which we believe to be worthwhile to be brought to the reader's attention. There has been a lot of papers where penalization of some goodness-of-fit criterion was proposed with rigorous proofs of consistency or some other guarantee of the resulting estimator of the true number of components. In [42] the penalized NPMLE was considered with a penalization function α n,m satisfying α n,m+1 > α n,m and lim sup n→∞ α n,m /n = 0. Under some regularity conditions on the mixture model, it is very rigorously shown that the estimator is at least equal to the true number of components as the sample n → ∞ with probability 1. The method requires only computation of the NPMLE for a number of values of m, which can be done using for example a support reduction algorithm as described in [70] or [33]. As it is not known whether this penalized NPMLE is actually consistent, [18] constructed a penalized minimum-distance estimator which could be thought as a precursor of the minimum distance estimators reviewed in Sect. 4. Two main differences are to be noted though: Firstly, [18] consider distances between distribution functions instead densities. Secondly, the penalization function takes the form of −c n m j=1 log π j where π j , j = 1, . . . , m are the mixing probabilities and (c n ) n a sequence converging to 0 as n → ∞. Consistency of the penalized minimum distance estimator of the true complexity is shown when one chooses c n such that the distance between the empirical distribution function and the true mixture distribution is O(c n ) almost surely. In the special case where one wants to decide between m 0 = 1 (homogeneity) and m 0 = 2, [17] propose a method based on modifying the likelihood ratio test. The modification operates first on the log-likelihood function by adding a negative penalty of the form C log(4π(1 − π)) with π > 0 the mixing probability and C > 0 some chosen constant. The penalty clearly discourages the MLE, under the null hypothesis, from fitting a mixing probability that is close to 0. The modification is motivated by the desire of overcoming the issues associated with the nesting F 1 ⊂ F 2 (boundary problem and non-uniqueness of representing the null hypothesis) already described in some details in Sect. 5. If π f φ 1 + (1 − π) f φ 2 is the mixture density, and if we denote by l n the modified log-likelihood, the the modified LRT is given by 2 l n (π,φ 1 ,π 2 ) − l n (1/2,φ,φ) , where (π,φ 1 ,φ 2 ) andφ are the MLE under the alternative and null hypotheses. One very interesting theoretical result is that the LRT is shown to converge weakly to mixture (1/2) : (1/2) of a Dirac at 0 and χ 2 (1) . This limit distribution can be then used to construct asymptotic critical region for rejecting homogeneity. As for the constant C, it is recommended to take log M if it is believed that φ 1 , φ 2 ∈ [−M, M] although the results do not seem to be too much sensitive to taking other values for C.
Bayesian approaches were also used to make inference about the number of components of a mixture. In this framework, this number is viewed as a random variable drawn from some prior distribution and the corresponding posterior distribution is then derived and subsequently used for inference purposes. For mixtures of Gaussian distributions whose means and variances are regarded as random variables drawn from a Dirichlet process, we refer to work of [25]. The posterior distribution can be approximated using Monte Carlo. In [60] the authors make use of reversible jump Markov chain Monte Carlo methods to conduct a a more general Bayesian analysis while restricting attention to Gaussian mixtures. In their analysis, the authors put themselves in the setting where no strong prior information on the components of such a mixture is available. In their work, the reader finds a thorough sensitivity analysis including the dependence of the posterior distribution of the number of components on the chosen prior for the means and variances. In the very interesting work of [16] the authors study the frequentist properties of Bayesian estimators of the true order in nested models including mixture models. Bounds on underestimation and overestimation of the Bayesian estimators are obtained. In particular, it is shown that, under some regularity conditions, the probability of underestimation decays exponentially. For further articles using Bayesian theory for clustering, we can refer to [51], [61], [30], [53] and the references therein.
We finish this section by drawing the reader's attention to existence of whole body of literature on mixture estimation and clustering, at the intersection of Statistics and Computing. This includes research papers on extensions or modifications of the famous EM-algorithm and numerical implementation of various information criteria. We refer to [15] where an entropy criteria was considered, which is derived from a simple relationship between the likelihood and classification likelihood of a mixture. In [27] an algorithm based on the Minimum Message Length (MML) criterion was implemented with the aim of selecting the best overall mixture model given the observed data (using a variant of the EM-algorithm). The very recent paper of [32] presents a novel form of cross-validation approach which is adaptive to the data. The paper contains an excellent literature review and the ideas discussed there are highly relevant, especially in connection with the question of how to choose the penalty function in the penalized methods reviewed above.

Some Conclusions
As acknowledged in the literature on the mixture model estimation and is supported by the results presented in Sect. 6, estimation of the number of components in a mixture distribution is a challenging task. None of the methods examined in the previous sections of the present survey can be regarded as a reliable universal tool that can be applied in any setting without second thoughts.
The widespread use among practitioners of BIC and ICL techniques is justified. These methods are easy to implement, do not require much computational time and outperform in many settings the other methods we have reviewed here. Whenever ICL tends to underestimate the number of components, which happens when they are not well separated, BIC does not seem to ehibit the same behavior producing more reliable estimates of the mixture order.
The LRT approach appears to be beneficial in terms of accuracy for small and average sample sizes, in particular in the settings when the mixture is not well-separated or when some components in the mixture noticeably dominate the other components or whenever the actual number of components in the mixture is large (e.g. 3 or 4 components). The disadvantage of this technique is the large amount of time needed for the computation.
The bootstrap methods on the average perform better than their penalty-based counterparts. In most of the settings MHD approach with bootstrap is more accurate than both MHD with the AIC-based penalty term and MHD with SBC-based term. The undeniable advantage of the procedures using the bootstrap is that the obtained estimates do not depend on the form of the penalty term, thus the errors resulting from the poor choice thereof can be avoided. The disadvantage however is the computational intensity of this procedure.
In the settings where the mixtures have well separated components LRT and MHD with bootstrap approaches provide an improvement over the other methods whenever the number of the components is more than 2. If a mixture comprises only 2 components and many observations are available, any of the methods can be applied.
The distance-based methods and the LRT approach can be used in the setting where the parameter values need to be estimated. If there is no such requirement, methods based on the Hankel matrix of moments of the mixing distributions can be used. For small sample sizes the scaled version of the Hankel matrix-based method achieves better performance than other methods. But generally it does not provide better accuracy than either the other methods, nor one can benefit significantly in terms of computational time.
Funding Open access funding provided by Swiss Federal Institute of Technology Zurich.

Conflict of interest statement
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.