Finite mixtures are a natural and ubiquitous class of models for modeling heterogeneous populations consisting of distinct subpopulations. Multiple volumes have been written regarding the application, analysis, and implementation of such models, including the volumes of Titterington et al. (1985); Lindsay (1995); Peel and McLachlan (2000a); Fruhwirth-Schnatter et al. (2019); Ng et al. (2019); Chen (2023), and Yao and Xiang (2024), among many others. When conducting mixture modeling, a frequently required task is that of order selection. In this text, we will take an information criteria (IC) approach following the tradition of Akaike (1974) as espoused in the texts of Anderson and Burnham (2004); Claeskens and Hjort (2008), and Konishi and Kitagawa (2008). We characterise the task of order selection and the IC method as follows.

Let \(\left( \varOmega ,\mathcal{A},\textrm{P}\right) \) be a probability space with expectation operator \(\textrm{E}\) and suppose that \(X:\varOmega \rightarrow \mathbb {X}\) is a random map with density function \(f_{0}:\mathbb {X}\rightarrow \mathbb {R}_{>0}\) with respect to dominating measure \(\mathfrak {m}\), in the sense that \(\textrm{P}\left( X\in \mathbb {A}\right) =\int _{\mathbb {A}}f_{0}\textrm{d}\mathfrak {m}\), for each \(\mathbb {A}\in \mathcal{B}\left( \mathbb {X}\right) \) (the Borel \(\sigma \)-algebra of \(\mathbb {X}\)). Let \(\phi :\mathbb {X}\times \mathbb {T}\rightarrow \mathbb {R}_{>0}\) be a probability density function (PDF) dependent on parameter \(\theta =\left( \vartheta _{1},\dots ,\vartheta _{m}\right) \in \mathbb {T}\subset \mathbb {R}^{m}\) in the sense that \(\int _{\mathbb {X}}\phi \left( x;\theta \right) \mathfrak {m}\left( \textrm{d}x\right) =1\), for each \(\theta \), and define the k-component mixture of \(\phi \)-densities as

where \(\left[ k\right] =\left\{ 1,\dots ,k\right\} \) and \(\psi _{k}=\left( \pi _{1},\dots ,\pi _{k},\theta _{1},\dots ,\theta _{k}\right) \in \mathbb {S}_{k}\subset \mathbb {R}^{\left( 1+m\right) k}\) (the parameter space of densities in \(\mathcal{M}_{k}^{\phi }\)). Here, \(\phi \) is referred to as the component PDF.

We suppose that \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\) but \(f_{0}\notin \mathcal{M}_{k_{0}-1}^{\phi }\) for some \(k_{0}\in \left[ \bar{k}\right] \), and some sufficiently large \(\bar{k}\in \mathbb {N}\). Here, the extra stipulation that \(f_{0}\notin \mathcal{M}_{k_{0}-1}^{\phi }\) is required because \(\left( \mathcal{M}_{k}^{\phi }\right) _{k\in \left[ \bar{k}\right] }\) is a nested sequence in the sense that \(\mathcal{M}_{k}^{\phi }\subset \mathcal{M}_{k+1}^{\phi }\), and thus if \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\), it also holds that \(f_{0}\in \mathcal{M}_{k_{0}+1}^{\phi }\), therefore the condition alone does not uniquely define \(k_{0}\), which we refer to as the true order of \(f_{0}\). We may interpret \(k_{0}\) as the order of the smallest (or most parsimonious) model \(\mathcal{M}_{k}^{\phi }\) that contains \(f_{0}\). Let \(\ell _{k}\left( \psi _{k}\right) =-\textrm{E}\left\{ \log f_{k}\left( X;\psi _{k}\right) \right\} \) denote the average negative log-density of \(f_{k}\left( \cdot ;\psi _{k}\right) \in \mathcal{M}_{k}^{\phi }\), for each \(k\in \left[ \bar{k}\right] \). Then, by virtue of the divergence property of the Kullback–Leibler divergence (cf. e.g., Csiszár (1995)), it also holds that

assuming that the minimum is well-defined in each case, which is a useful characterisation that we will exploit in the sequel.

Suppose that we observe a data set of n independent and identically distributed (IID) replicates of X: \(\textbf{X}_{n}=\left( X_{1},\dots ,X_{n}\right) \subset \mathbb {X}\). Then, the problem of order selection can be characterised as the construction of an estimator \(\hat{k}_{n}=\hat{k}\left( \textbf{X}_{n}\right) \) of \(k_{0}\) with desirable properties. In particular, we want \(\hat{k}_{n}\) to be consistent in the sense that

If we let \(\ell _{k,n}\left( \psi _{k}\right) =-n^{-1}\sum _{i=1}^{n}\log f_{k}\left( X_{i};\psi _{k}\right) \) denote the empirical average negative log-likelihood of \(f_{k}\left( \cdot ;\psi _{k}\right) \), the natural estimator of \(\ell _{k}\left( \psi _{k}\right) \), then the method of IC suggests that we construct the estimator of \(k_{0}\) in the form

for some appropriately chosen penalties \(\textrm{pen}_{k,n}\). Here, we may recognise \(\min _{\psi _{k}\in \mathbb {S}_{k}}\ell _{k,n}\left( \psi _{k}\right) \) as the negative maximum likelihood for the mixture model of order k.

We can easily identify this approach with some traditional methods, such as the Akaike information criterion (AIC; Akaike (1974)) and the Bayesian information criterion (BIC; Schwarz (1978)), by taking \(\textrm{pen}_{k,n}^{\textrm{AIC}}=d_{k}n^{-1}\) and \(\textrm{pen}_{k,n}^{\textrm{BIC}}=d_{k}\left( 2n\right) ^{-1}\log n\), respectively, where \(d_{k}\) denotes the parameter count used in the criterion. Let \(\textrm{Ln}\left( n\right) =\log \left( e\vee n\right) \) denote the logarithm function truncated to 1, from below, where \(a\vee b=\max \left\{ a,b\right\} \), and denote its \(\nu \)-fold composition by \(\textrm{Ln}^{\circ \nu }\left( n\right) \) (e.g., \(\textrm{Ln}^{\circ 3}\left( n\right) =\textrm{Ln}\circ \textrm{Ln}\circ \textrm{Ln}\left( n\right) \)). In this work, we introduce the \(\nu \)-BIC and \(\epsilon \)-BIC families of penalties of the forms \(\textrm{pen}_{k,n}^{\nu }=\alpha \left( k\right) n^{-1}\textrm{Ln}^{\circ \nu }\left( n\right) \log n\) and \(\textrm{pen}_{k,n}^{\epsilon }=\alpha \left( k\right) n^{-1}\left( \log n\right) ^{1+\epsilon }\), for any choice of \(\epsilon >0\) and \(\nu \in \mathbb {N}\), and for \(\alpha :\mathbb {N}\rightarrow \mathbb {R}_{>0}\) strictly increasing. A BIC-style calibration is obtained by taking \(\alpha \left( k\right) =d_{k}/2\). We shall demonstrate that these modest modifications of the BIC are consistent under weaker assumptions than those currently available for BIC consistency in mixtures, as we shall elaborate upon in the sequel. We also provide a complementary misspecification result: when \(f_{0}\) need not belong to any \(\mathcal{M}_{k}^{\phi }\), any IC satisfying a mild vanishing-penalty condition will eventually select an order whose best-fitting mixture model is Kullback–Leibler optimal among the candidate orders.

For regular models, order selection consistency of the BIC was proveds in the works of Nishii (1988); Vuong (1989), and Sin and White (1996). Unfortunately, the mixture model setting is not regular and thus alternative proofs of consistency are required. A prime and often cited example of such a result is that of Keribin (2000) who establish the consistency of the BIC for mixture models under strong assumptions. Namely, among other things, the proof requires that the component density have up to five derivatives and that all of these higher derivatives have finite third moments. A remarkable result of Gassiat and Van Handel (2012) (see also Gassiat 2018, Sec. 4.3) shows that the BIC is order consistent for choosing between models when \(\bar{k}=\infty \), albeit at the price of much stronger assumptions. In particular, even when restricted to the class of location mixtures, the proof requires that \(\phi \) has up to three derivatives and that all derivatives up to second order have finite third moments, while third order derivatives have second moments.

Beyond the BIC, Drton and Plummer (2017) proposed the so-called singular BIC (sBIC) as an approach that operationalises the expansion \(n\ell _{k,n}\left( \psi _{k}\right) +\lambda _{k}\left( f_{0}\right) \log n+\left\{ m_{k}\left( f_{0}\right) -1\right\} \log \log n+O_{\textrm{P}}\left( 1\right) \) of the so-called Bayes free energy of Watanabe (2013), related to the widely-applicable BIC (WBIC), where the coefficients \(\left( \lambda _{k}\left( f_{0}\right) ,m_{k}\left( f_{0}\right) \right) _{k\in \left[ \bar{k}\right] }\) are dependent on the true and unknown model \(f_{0}\). As noted in Chen (2023, Sec. 16.3), the specification of the sBIC arises as a solution to a nonlinear system of stochastic equations making it challenging to faithfully implement in practice. In Baudry (2015), the log-likelihood is replaced by the so-called complete data log-likelihood and a consistency result is obtained with the usual BIC-form penalty. However, the proof makes the requirement that the limiting objective functionals of all order \(k\in \left[ \bar{k}\right] \) have positive definite Hessian, which is not possible in the context of mixture models, as typically if \(k_{0}<k\), then the minimisers of any risk over \(\mathcal{M}_{k}^{\phi }\) will have connected components. In the context of test-based order selection, this phenomenon is discussed in Chen (2023, Ch. 9).

Sacrificing the requirement that \(\textrm{pen}_{k,n}=\tilde{O}\left( n^{-1}\right) \) (where we use the soft-Oh notation \(\tilde{O}\) to denote the Landau order up to logarithmic factors; cf. Cormen et al. (2002), Sec. 3-5), Nguyen (2024) propose to take penalties such that \(\lim _{n\rightarrow \infty }\sqrt{n}\textrm{pen}_{k,n}=\infty \), e.g. \(\textrm{pen}_{k,n}=\tilde{O}\left( n^{-1/2}\right) \), within their PanIC framework for model selection in general risk minimisation problems (see also Westerhout et al. (2024)). Using these larger penalties, it is possible to prove model selection consistency under very mild conditions, with data potentially dependent, and where the objective functions (including but not necessarily log-likelihoods) may be non-differentiable or even discontinuous and unmeasurable. Although these results are very general, it is clear that for sufficiently well-behaved models, \(\tilde{O}\left( n^{-1/2}\right) \) penalties are too large and as a consequence, the corresponding IC \(\hat{k}_{n}\) may be inefficient and thus larger than necessary values of n are required for consistency to manifest, in practice.

In this work, we will prove that the slightly larger penalties \(\textrm{pen}_{k,n}^{\nu }\) and \(\textrm{pen}_{k,n}^{\epsilon }\) allow for a substantial relaxation of the conditions required for order consistency. In particular, we will show that consistency can be proved when the parameter space \(\mathbb {S}_{k}\) is compact, the log-density class is a Glivenko–Cantelli (GC) class, and the component density functions are Lipschitz over \(\mathbb {T}\) with an \(\mathcal {L}_{1}(\mathfrak m)\) envelope. These conditions are far less onerous than those of Keribin (2000) and allow for application in cases where the component densities \(\phi \) need not even be differentiable. In fact, penalties \(\textrm{pen}_{k,n}^{\epsilon }\) were considered within the framework of Keribin (2000) but the larger penalisation was not exploited to weaken the regularity conditions.

In fact, for \(\nu \) chosen sufficiently large or \(\epsilon \) chosen sufficiently small, the decision made by using the \(\nu \)-BIC or \(\epsilon \)-BIC versus that by the BIC will be the same, in practical settings. For instance, when \(\nu =3\), \(\textrm{Ln}^{\circ \nu }\left( n\right) \le 1\) for all \(n\le \exp ^{\circ 3}\left( 1\right) \approx 3.8\times 10^{6}\) and \(\textrm{Ln}^{\circ \nu }\left( n\right) \le 1.1\) for all \(n\le \exp ^{\circ 3}\left( 1.1\right) \approx 5.7\times 10^{8}\). Similarly, taking \(\epsilon =0.02\), \(\left( \log n\right) ^{0.02}\le 1.1\) for all \(n\le \exp \left( 117.3909\right) \approx 9.6\times 10^{50}\). Thus, one key interpretation of our results is that they provide theoretical support for penalties that are numerically very close to the BIC outside the regularity framework of Keribin (2000). In finite sample studies, the use of the BIC may therefore be viewed as numerically close to the use of \(\nu \)-BIC or \(\epsilon \)-BIC penalties, up to a negligible multiplicative factor. We also use our framework to clarify two limitations of consistent IC-based order selection in mixtures. First, within the general sufficient conditions used to establish consistency, there is no universally minimal BIC-scale penalty, so consistency alone does not single out an optimal IC. Second, we establish a tension between order consistency and minimax optimality in mixtures in the style of Yang (2005). In particular, we show that in a simple Gaussian mixture family, an AIC-like criterion achieves the parametric minimax Hellinger risk, whereas any order-consistent criterion (including the BIC and our proposals) must be minimax-inefficient by an unbounded factor. To the best of our knowledge, this mixture model AIC/BIC tension in mixture models has not previously been stated explicitly.

Indeed, the method of IC is not the only approach for order selection and estimation under model uncertainty in mixture models. For example, model selection using hypothesis testing has been considered in McLachlan (1987); Polymenis and Titterington (1998); Wasserman et al. (2020); Nguyen et al. (2023), and the many approaches covered in Chen (2023). Shrinkage estimator approaches have been proposed in Chen and Khalili (2009); Huang et al. (2017, 2022), and Budanova (2025). And bounds of risk functions for estimation under model uncertainty have been considered by Li and Barron (1999); Rakhlin et al. (2005); Klemelä (2007); Maugis and Michel (2011); Manole and Ho (2022), and Chiu Chong et al. (2024), among many others. We do not comment or compare these disparate approaches to model selection, and direct the reader to dedicated texts, such as Peel and McLachlan (2000a, Ch. 6), McLachlan and Rathnayake (2014), Fruhwirth-Schnatter et al. (2019, Ch. 7), and Chen (2023, Ch. 16).

To conclude, we note that beyond mixtures of density functions, one can consider IC for conditional density models, as well. Results in this setting include the works of Naik et al. (2007); Hafidi and Mkhadri (2010); Depraetere and Vandebroek (2014), and Hui et al. (2015). It is possible to adapt our approach to provide inference in this setting, and we will consider such situations as an example application of our theory.

The remainder of the manuscript is organised as follows. Section 2 introduces the notational conventions used throughout the paper, presents the technical preliminaries, and states our main consistency results, including a misspecification result in which order selection targets Kullback–Leibler optimal approximating models. Section 3 provides illustrative examples that demonstrate applications of our theory. Section 4 then discusses limitations of consistent model selection in mixture models, including penalty flexibility and the aforementioned AIC/BIC tension. Concluding remarks and discussion are then made in Sect. 5. Complete proofs of our technical results and primary theorems are provided in the supplementary material. Auxiliary assumptions used for comparison are provided in the Appendix.

We retain the setup and notation from Sect. 1 and follow the usual notational conventions for empirical processes (see, e.g., van de Geer (2000)). In particular, \(\mathcal{M}_{k}^{\phi }\) denotes the class of \(k\)-component mixture densities, \(f_{k}(\cdot ;\psi _{k})\) denotes a density in that class, \(\mathbb {S}_{k}\) denotes its parameter space, \(\ell _{k}\) and \(\ell _{k,n}\) denote the population and empirical average negative log-likelihoods, \(\textrm{pen}_{k,n}\) denotes the order-\(k\) penalty, and \(\hat{k}_{n}\) denotes the resulting IC-selected order. Let \(P\) denote the true probability measure of \(X\), and use empirical-process notation \(Pf=\textrm{E}f(X)\) for an integrable function \(f\). Let \(P_n\) denote the empirical measure, so that \(P_nf=n^{-1}\sum _{i=1}^n f(X_i)\). For a functional class \(\mathcal{F}\), \(\ell ^\infty (\mathcal{F})\) is the space of bounded functionals on \(\mathcal{F}\) equipped with the uniform norm \(\Vert h\Vert _{\mathcal{F}}=\sup _{f\in \mathcal{F}}|h(f)|\). When the class is clear from context, \(P_n-P\) is viewed as an element of \(\ell ^\infty (\mathcal{F})\). A class \(\mathcal{F}\) is GC if it satisfies the uniform strong law of large numbers \(\Vert P_n-P\Vert _\mathcal {F}\xrightarrow [n\rightarrow \infty ]{\mathrm {a.s.}}0\). We write \(a_n\asymp b_n\) if there exist constants \(0<c\le C<\infty \) and \(n_0\in \mathbb {N}\) such that \(cb_n\le a_n\le Cb_n\) for all \(n\ge n_0\), and write \(g(n)\ll h(n)\) if \(g(n)/h(n)\rightarrow 0\).

For \(p\ge 1\), \(\mathcal{L}_p(P)=\mathcal{L}_p(\mathbb {X},\mathcal{B}(\mathbb {X}),P)\) denotes the usual Lebesgue space, with its associated norm. For a vector \(\psi \in \mathbb {R}^q\), \(\Vert \psi \Vert \) denotes the Euclidean norm and \(\Vert \psi \Vert _p\) denotes the \(\ell _p\) norm. For densities \(f\) and \(g\) with respect to \(\mathfrak {m}\), the Hellinger distance is

For probability measures \(Q\) and \(R\) on \((\mathbb {X},\mathcal{B}(\mathbb {X}))\), the Kullback–Leibler divergence is

If \(Q\) and \(R\) have densities \(f\) and \(g\) with respect to \(\mathfrak {m}\), we also write

We write

for the corresponding minimized risks. Thus \(L_k\) and \(\hat{L}_{k,n}\) denote the minimized population and empirical risks, respectively, as distinct from the parameter-indexed risks \(\ell _k\) and \(\ell _{k,n}\). When likelihood comparisons are written in maximisation form over density classes \(\mathcal {F}_k\), we use, for any probability measure \(Q\) for which the displayed quantity is well defined,

When the relevant law is clear, we write \(\varLambda _k\) in place of \(\varLambda _k(Q)\). Consequently, when \(\mathcal {F}_k=\mathcal{M}_k^\phi \), one has \(\hat{L}_{k,n}=-\varLambda _{k,n}\) and \(L_k=-\varLambda _k(P)\).

For a norm \(\Vert \cdot \Vert \), a collection \(( [f_j^{\mathrm L},f_j^{\mathrm U}] )_{j\in [N]}\) is a \(\delta \)-bracketing of \(\mathcal{F}\) if, for each \(f\in \mathcal{F}\), there is a \(j\) such that \(f_j^{\mathrm L}\le f\le f_j^{\mathrm U}\) and \(\Vert f_j^{\mathrm U}-f_j^{\mathrm L}\Vert \le \delta \), for every \(j\in [N]\). The bracketing number \(N_{\left[ \right] }(\delta ,\mathcal{F},\Vert \cdot \Vert )\) is the smallest such \(N\), with value \(\infty \) if no finite bracketing exists, and \(H_{\left[ \right] }(\delta ,\mathcal{F},\Vert \cdot \Vert )=\log N_{\left[ \right] }(\delta ,\mathcal{F},\Vert \cdot \Vert )\) is the corresponding bracketing entropy. When \(\mathcal{F}\subset \mathcal{L}_p(P)\) and \(\Vert \cdot \Vert \) is the \(\mathcal{L}_p(P)\) norm, we write \(N_{\left[ \right] }(\delta ,\mathcal{F},\mathcal{L}_p(P))\) for short. For a metric space \((\mathcal{G},d)\), a collection \((g_j)_{j\in [N]}\) is a \(\delta \)-covering of \(\mathcal{G}\) if each \(g\in \mathcal{G}\) is within distance \(\delta \) of some \(g_j\); the smallest such \(N\) is denoted \(N(\delta ,\mathcal{G},d)\).

With this notation, the following results are useful for proving our technical outcomes.

If \(\textbf{X}_{n}\) is IID and \(\mathcal{F}\subset \mathcal{L}_{1}\left( P\right) \) is a class of measurable functions admitting an envelope \(F\in \mathcal{L}_{1}\left( P\right) \) and \(N_{\left[ \right] }\left( \delta ,\mathcal{F},\mathcal{L}_{1}\left( P\right) \right) <\infty \), for every \(\delta >0\), then \(\mathcal{F}\) is a GC class. In particular, if

\(\left( \mathcal{G},d\right) \) is compact and \(g\mapsto f_{g}\left( x\right) \) is continuous for P-almost all \(x\in \mathbb {X}\), then \(N_{\left[ \right] }\left( \delta ,\mathcal{F},\mathcal{L}_{1}\left( P\right) \right) <\infty \), for every \(\delta >0\), if there exists an envelope function \(F\in \mathcal{L}_{1}\left( P\right) \), such that \(\sup _{g\in \mathcal{G}}\left| f_{g}\right| \le F\).

See van de Geer (2000, Thm. 2.4.1 and Lem. 3.10). \(\square \)

As an alternative to the bracketing entropy condition, above, one can also verify that \(\mathcal{F}\) is a GC class via random covering entropy methods (see, e.g., van de Geer (2000), Sec. 3.6). We do not use such conditions in this text since bracketing entropy bounds pair better with the rest of our theoretical approach as shall be seen, in the sequel.

Throughout this text, we shall assume that all random quantities of interest are measurable. Indeed, all of the classes \(\mathcal{F}\) considered in this work can be understood as classes of Carathéodory integrands that are indexed by Euclidean spaces and thus are measurable, with measurable maxima/minima, sums, compositions, and other required manipulations (cf. Rockafellar and Wets (2009), Ch. 14).

If \(\mathcal{F}\) is parameterised as per (1), such that there exists an envelope function \(F:\mathbb {X}\rightarrow \mathbb {R}_{>0}\) satisfying

then \(N_{\left[ \right] }\left( 2\delta \left\| F\right\| ,\mathcal{F},\left\| \cdot \right\| \right) \le N\left( \delta ,\mathcal{G},d\right) \). Furthermore, if \(\mathcal{G}=\mathbb {S}\subset \mathbb {R}^{q}\) is compact and d is the Euclidean distance, then there exists a constant \(K>0\) such that

for every \(0<\delta <\textrm{diam}\left( \mathbb {S}\right) \).

See van der Vaart and Wellner (van der Vaart and Wellner (2023), Thm. 2.7.17) and van der Vaart (2000, Example 19.7). \(\square \)

Next, suppose that \(\mathcal{F}\) is a subset of probability density functions on \(\mathbb {X}\), in the sense that for each \(f\in \mathcal{F}\), \(f\left( x\right) >0\) for each \(x\in \mathbb {X}\) and \(\int _{\mathbb {X}}f\textrm{d}\mathfrak {m}=1\). For each \(n\in \mathbb {N}\), we let \(\hat{f}_{n}:\varOmega \rightarrow \mathcal{F}\) be a maximum likelihood estimator (MLE), defined by

where the existence of MLEs is guaranteed in our application. Using the Hellinger distance defined in Sect. 2.1, write the subset of \(\mathcal{F}\) that lies within the \(\delta \)-Hellinger ball of \(f_{0}\), for \(0<\delta \le 1\), as \(\mathcal{F}\left( \delta \right) =\left\{ f\in \mathcal{F}:\mathfrak {h}\left( f,f_{0}\right) \le \delta \right\} \). Further define the following entropy integral:

The next result, our main technical tool, is adapted from van de Geer (2000, Cor. 7.5) (see the supplementary material).

If \(f_{0}\in \mathcal{F}\), \(\textbf{X}_{n}\) is IID and there is some \(\varPsi :\mathbb {R}_{>0}\rightarrow \mathbb {R}_{>0}\) such that \(\varPsi \left( \delta \right) \ge J\left( \delta \right) \), where \(\varPsi \left( \delta \right) /\delta ^{2}\) is a non-increasing function of \(\delta \), then, for some constant \(c_{1}>0\),

for every \(\delta \ge \delta _{n}\), such that \(\delta _{n}\) satisfies \(\sqrt{n}\delta _{n}^{2}\ge c\varPsi \left( \delta _{n}\right) \), for a universal constant \(c>0\).

We state the following assumptions, repeating some previously stated items, for convenience.

A1:

\(X\in \mathbb {X}\) has PDF \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\setminus \mathcal{M}_{k_{0}-1}^{\phi }\), for some \(k_{0}\in \left[ \bar{k}\right] \) (with the convention \(\mathcal{M}_{0}^{\phi }=\varnothing \)), and \(\textbf{X}_{n}=\left( X_{1},\dots ,X_{n}\right) \) is a sample of IID replicates of X.

A2:

\(\phi \) is strictly positive; i.e., \(\phi \left( x;\theta \right) >0\), for each \(x\in \mathbb {X}\) and \(\theta \in \mathbb {T}\subset \mathbb {R}^{m}\), where \(\mathbb {T}\) is compact.

A3:

\(\phi \) is a Carathéodory integrand in the sense that \(\phi \left( x;\cdot \right) :\mathbb {T}\rightarrow \mathbb {R}\) is continuous for each fixed \(x\in \mathbb {X}\), and \(\phi \left( \cdot ;\theta \right) :\mathbb {X}\rightarrow \mathbb {R}\) is measurable, for each fixed \(\theta \in \mathbb {T}\).

A4:

There exists a \(G_{1}\in \mathcal{L}_{1}\left( P\right) \) such that, for every \(x\in \mathbb {X}\),

$$ \max _{\theta \in \mathbb {T}}\left| \log \phi \left( x;\theta \right) \right| <G_{1}\left( x\right) . $$

A5:

There exist functions \(L_{\phi }:\mathbb {X}\rightarrow \mathbb {R}_{\ge 0}\) and \(G_{2}\in \mathcal{L}_{1}\left( \mathfrak {m}\right) \) such that, for every \(\theta _{1},\theta _{2}\in \mathbb {T}\) and \(x\in \mathbb {X}\),

$$\begin{aligned} \left| \phi \left( x;\theta _{1}\right) -\phi \left( x;\theta _{2}\right) \right| \le L_{\phi }\left( x\right) \left\| \theta _{1}-\theta _{2}\right\| _{1}, \end{aligned}$$

(4)

where \(\max _{\theta \in \mathbb {T}}\phi \left( x;\theta \right) +L_{\phi }\left( x\right) \le G_{2}\left( x\right) \).

When \(\mathbb {T}\) is convex and \(\phi \left( x,\cdot \right) :\mathbb {T}\rightarrow \mathbb {R}\) is continuously differentiable, for each \(x\in \mathbb {X}\), we can replace A5 with the following condition.

A5*:

The set \(\mathbb {T}\) is convex, and there exists a function \(G_{2}\in \mathcal{L}_{1}\left( \mathfrak {m}\right) \) such that, for each \(x\in \mathbb {X}\),

$$ \max _{\theta \in \mathbb {T}}\phi \left( x;\theta \right) +\max _{j\in \left[ m\right] }\max _{\theta \in \mathbb {T}}\left| \frac{\partial \phi \left( x;\theta \right) }{\partial \vartheta _{j}}\right| \le G_{2}\left( x\right) . $$

We also impose the following penalty conditions, under which penalties are non-negative throughout.

B1:

For each \(k\in \left[ \bar{k}\right] \), \(\textrm{pen}_{k,n}\ge 0\) for all n and \(\lim _{n\rightarrow \infty }\textrm{pen}_{k,n}=0\).

B2:

For each \(1\le k<l\le \bar{k}\),

$$ \lim _{n\rightarrow \infty }\frac{n}{\log n}\left\{ \textrm{pen}_{l,n}-\textrm{pen}_{k,n}\right\} =\infty . $$

Let us provide some reasoning for the assumptions. A1 identifies that the true PDF \(f_{0}\) is within the class of mixture models \(\mathcal{M}_{k}^{\phi }\), where \(k\le \bar{k}\). We also have the availability of an IID sample of n replicates \(\textbf{X}_{n}\) from the true probability model, required for application of Lemma 1 and Proposition 1. A2 provides the necessary compactness of the parameter space for bounding the metric entropy using Lemma 2 and ensures that the log-densities \(\log \phi \left( x;\theta \right) \ne \infty \), for each pair of x and \(\theta \). A3 implies that the density expressions, log-density expressions, and their minima/maxima and averages remain measurable. A4 provides envelope functions so that we can apply Lemma 1 to prove the convergence of the average log-likelihoods to their limits, for each k. And A5 provides the required envelope functions for bounding the bracketing entropy using Lemma 2. Note that when \(\phi \) is continuously differentiable in \(\theta \), for each fixed x, and \(\mathbb {T}\) is convex, A5* implies A5 by the mean value theorem. B1 imposes non-negative vanishing penalties and guarantees that the penalty does not interfere with comparisons between models \(\mathcal{M}_{k}^{\phi }\) and \(\mathcal{M}_{k_{0}}^{\phi }\), with respect to their maximum average log-likelihoods, when \(k<k_{0}\). B2 then provides sufficient penalisation to distinguish between \(\mathcal{M}_{k}^{\phi }\) and \(\mathcal{M}_{k_{0}}^{\phi }\), with respect to their model complexities, when \(k>k_{0}\).

The following pair of lemmas, proved in the supplementary material, constitute the main technical contribution of the text and will imply the main result.

Under A1–A4, for each \(k\in \left[ \bar{k}\right] \),

Under A1–A3 and A5, for each \(k,l\ge k_{0}\),

As a consequence, we have the following generic consistency result and our main result as a corollary.

Under A1–A5, B1, and B2, \(\lim _{n\rightarrow \infty }\textrm{P}\left( \hat{k}_{n}=k_{0}\right) =1\).

Write \(\hat{L}_{k,n}=\min _{\psi _{k}\in \mathbb {S}_{k}}\ell _{k,n}(\psi _{k})\) and \(L_{k}=\min _{\psi _{k}\in \mathbb {S}_{k}}\ell _{k}(\psi _{k})\). Since \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\setminus \mathcal{M}_{k_{0}-1}^{\phi }\) and the model sequence is nested, \(L_{k}>L_{k_{0}}\) for every \(k<k_{0}\). Indeed, \(L_{k}-L_{k_{0}}=\inf _{f\in \mathcal{M}_{k}^{\phi }}\mathfrak {K}(P\Vert P_{f})\), where \(P_{f}\) denotes the measure with density f. If the infimum were zero, Pinsker’s inequality and compactness of the parameter space, together with the \(\mathcal{L}_{1}(\mathfrak {m})\)-continuity supplied by A5, would imply that \(f_{0}\in \mathcal{M}_{k}^{\phi }\), contradicting A1. Because there are finitely many such k, there exists \(\eta >0\) such that \(L_{k}-L_{k_{0}}>2\eta \) for all \(k<k_{0}\). Lemma 3 and B1 imply that, with probability tending to one, \(\hat{L}_{k,n}+\textrm{pen}_{k,n}>\hat{L}_{k_{0},n}+\textrm{pen}_{k_{0},n}\) for every \(k<k_{0}\).

For \(k>k_{0}\), Lemma 4 gives

whereas B2 gives

Hence, with probability tending to one, \(\hat{L}_{k,n}+\textrm{pen}_{k,n}>\hat{L}_{k_{0},n}+\textrm{pen}_{k_{0},n}\) for every \(k>k_{0}\). Since \(\bar{k}\) is finite, combining the underfitting and overfitting comparisons yields \(\textrm{P}(\hat{k}_{n}=k_{0})\rightarrow 1\). \(\square \)

Under A1–A5, for any strictly increasing \(\alpha :\mathbb {N}\rightarrow \mathbb {R}_{>0}\), \(\nu \in \mathbb {N}\) and \(\epsilon >0\), the \(\nu \)-BIC

and \(\epsilon \)-BIC

are consistent estimators of \(k_{0}\).

For every fixed \(\nu \in \mathbb {N}\) and \(\epsilon >0\), both \(\textrm{Ln}^{\circ \nu }\left( n\right) \log n=o(n)\) and \(\log ^{1+\epsilon }n=o(n)\). Thus \(\textrm{pen}_{k,n}^{\nu }\) and \(\textrm{pen}_{k,n}^{\epsilon }\) satisfy B1. For \(k<l\), observe that

and

as required, since \(\log \) and positive powers are increasing, and since \(\alpha \) is strictly increasing. \(\square \)

To illustrate the advantage of Theorem 1 over the main result of Keribin (2000) (i.e., Theorem 2.1), it is helpful to compare the assumptions, which we supply for convenience in the Appendix. Firstly, Keribin (2000) shares A1–A3, and makes the additional condition (Id) requiring identifiability regarding mappings between the functional spaces \(\mathcal{M}_{k}^{\phi }\) and the parameter spaces \(\mathbb {S}_{k}\), which is unnecessary in our approach. Condition (P1-a) then requires that there exists a \(G\in \mathcal{L}_{1}\left( P\right) \) such that \(\left| \log f\right| \le G\), for every \(f\in \mathcal{M}_{\bar{k}}^{\phi }\), which is implied by A4. Condition (P1-b) then requires that the component density \(\phi \) possesses partial derivatives up to order 5, where all the partial derivatives of orders 2, 3, and 5, are each enveloped over \(\mathbb {T}\) by integrable functions in \(\mathcal{L}_{3}\left( \mathfrak {m}\right) \). This compares directly with A5, where we assume that \(\phi \) need not even have one derivative, only requiring that it is Lipschitz, with Lipschitz constant in \(\mathcal{L}_{1}\left( \mathfrak {m}\right) \) and that \(\phi \) itself is enveloped by an \(\mathcal{L}_{1}\left( \mathfrak {m}\right) \) function. This substantially simplifies the assumptions and makes our approach useful in many situations where the method of Keribin (2000) does not apply.

Conditions (P2) and (P3), which our theory does not require, then make linear independence assumptions between the first and second order partial derivatives of \(\phi \), and require that a class of linear combinations of zeroth- and first-order derivatives of \(\phi \) be a Donsker class of functionals with continuous sample paths, respectively. Finally, in the average-risk normalisation used here, Condition (C1) requires positive, increasing, vanishing penalties and imposes the true-order separation condition

for each \(k_{0}<k\le \bar{k}\). Thus B2 strengthens the separation scale in Condition (C1) by a logarithmic factor. We believe this is a small price to pay for the reduction or elimination of the other conditions. In this setting we only require an infinitesimally larger penalty than that of the BIC, which is negligible in practice.

A complementary comparison can be made to the PanIC framework of Nguyen (2024), which treats model selection problems in a general empirical-risk minimisation form. For convenience, we reproduce in the Appendix the assumptions (N-A1)–(N-A2) and the penalty conditions (N-B1)–(N-B2) of Nguyen (2024). Specialising their setting to likelihood-based order selection by identifying \(\mathbb {T}_{k}\) with \(\mathbb {S}_{k}\) and taking the loss to be \(\ell _{k}(x;\psi _{k})=-\log f_{k}(x;\psi _{k})\), the hypotheses (N-A1)–(N-A2) amount to compactness of \(\mathbb {S}_{k}\) together with a global Lipschitz condition on the negative log-likelihood with an \(\mathcal{L}_{2}(P)\) envelope, which are slightly simpler to verify than A2–A5.

Under (N-A1)–(N-A2), Nguyen (2024, Thm. 1) establishes consistency for any penalised estimator whose penalties satisfy (N-B1)–(N-B2), where (N-B2) requires penalty size condition \(\sqrt{n}\{\textrm{pen}_{l,n}-\textrm{pen}_{k,n}\}\rightarrow \infty \) for \(k<l\), implying that \(\textrm{pen}_{l,n}-\textrm{pen}_{k,n}\) must dominate \(O\left( n^{-1/2}\right) \). By contrast, in the finite mixture likelihood setting, Lemma 4 shows that the relevant overfitting fluctuations are of order \(O\left( \log n/n\right) \). Consequently, the separation condition B2 in Theorem 1 only requires \(\textrm{pen}_{l,n}-\textrm{pen}_{k,n}\) to dominate \(O\left( \log n/n\right) \), permitting penalties that are asymptotically much smaller than those required by (N-B2) while delivering the same consistency guarantee that \(\textrm{P}(\hat{k}_{n}=k_{0})\rightarrow 1\) under A1–A5.

As discussed in Sect. 1, the penalty induced by the \(\epsilon \)-BIC and \(\nu \)-BIC can be made indistinguishable from that of the BIC when \(\epsilon \) is small or when \(\nu \) is large. In Table 1, we provide a schedule for the ratio of the rates of the \(\epsilon /\nu \)-BIC to that of the BIC, computed as \(\textrm{pen}_{k,n}^{\epsilon }/\textrm{pen}_{k,n}^{\textrm{BIC}}\) or \(\textrm{pen}_{k,n}^{\nu }/\textrm{pen}_{k,n}^{\textrm{BIC}}\), for the same choice of \(\alpha \), where \(\textrm{pen}_{k,n}^{\textrm{BIC}}=\alpha \left( k\right) n^{-1}\log n\). We observe that for sufficiently small choices of \(\epsilon \) or large choices of \(\nu \), the penalty ratios can be made to stay close to one for a larger range of sample sizes n.

It would appear that to improve upon the fluctuation bound \(O\left( \log n/n\right) \) in Lemma 4 (equivalently, the normalising rate \(n/\log n\), and subsequently B2), sharper bounds on the entropy \(H_{\left[ \right] }\left( u,\bar{\mathcal{F}}^{1/2}\left( \delta \right) ,\mathcal{L}_{2}\left( \mathfrak {m}\right) \right) \) are required that make use of the localising condition \(\mathfrak {h}\left( f,f_{0}\right) \le \delta \) of \(\bar{\mathcal{F}}^{1/2}\left( \delta \right) \). It is in the effort of obtaining such sharp bounds that necessitates the additional assumptions regarding the partial derivatives of \(\phi \) employed by Keribin (2000) and Gassiat and Van Handel (2012). Indeed, van de Geer (2000, Sec. 7.5) notes that in the finite dimensional parametric setting, one cannot obtain \(O\left( n\right) \) rates without exploiting the localising condition. On the other hand, in Theorem 1, the penalty need only satisfy B1 and B2, and thus need only go to zero and scale at a rate larger than \(O\left( n^{-1}\log n\right) \). In fact the two choices that define the \(\nu \)-BIC and \(\epsilon \)-BIC are arbitrary and were chosen because they can be made arbitrarily close to the BIC in practice. However, one may choose any sufficiently large penalty. For example, in Nguyen (2024) the choice of \(\textrm{pen}_{k,n}=O\left( n^{-1/2}\sqrt{\textrm{Ln}^{\circ \nu }\left( n\right) }\right) \) is proposed, for \(\nu \ge 1\), in conforming with the suggestion in the treatment of Sin and White (1996), and applied successfully to mixture model settings in Nguyen (2024) and Westerhout et al. (2024), within the PanIC framework.

In Corollary 1, the choice of \(\alpha \) need only be a strictly increasing function in k and as such any such choice is valid for the conclusion to hold. For example, one may ignore the multiplier due to m and simply take \(\alpha \left( k\right) =k\), or one may take \(\alpha \left( k\right) \) to be proportional to the number of free parameters in \(\mathcal{M}_{k}^{\phi }\), as in the usual BIC. Concretely, since \(\psi _{k}\) contains \((k-1)\) free mixing weights and km component parameters, a natural choice is

or the simplified version \(\alpha \left( k\right) =(m+1)k/2\) which we propose throughout, which ignores the simplex constraint. This indifference to the count of the number of so-called free parameters shows that one need not be so judicious regarding what defines a parameter, as in Fraley and Raftery (2002) and Baudry et al. (2010). Other admissible choices include \(\alpha \left( k\right) =k\) (counting mixture weights only; as considered in Keribin (2000)), \(\alpha \left( k\right) =km/2\) (counting component parameters only), or more aggressively \(\alpha \left( k\right) =k\log k\) or \(\alpha \left( k\right) =k^{\gamma }\) for \(\gamma >1\). Such freedom mirrors that in other consistent IC frameworks; for instance, the deterministic SWIC construction of Nguyen (2024, Rem. 3) permits any strictly increasing positive sequence to weight the common rate. More generally, the PanIC definition allows random penalty arrays, provided the corresponding stochastic penalty conditions are satisfied.

The results of Sect. 2.4 depend crucially on Assumption A1, which requires that \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\) for some \(k_{0}\in \left[ \bar{k}\right] \). In practice, however, it is not always known that \(f_{0}\) belongs to the family of mixtures under study (and in particular, a “true order” \(k_{0}\) need not exist). In such a case, the conclusions of Theorem 1 and Corollary 1 are no longer meaningful as statements about recovering \(k_{0}\), but the information criteria considered above still admit a sensible population interpretation. We therefore replace A1 by the following assumption.

A1*:

\(X\in \mathbb {X}\) has PDF \(f_{0}\) with respect to \(\mathfrak {m}\), and \(\textbf{X}_{n}=\left( X_{1},\dots ,X_{n}\right) \) is a sample of IID replicates of X.

Recall the minimized risks \(L_{k}\) and \(\hat{L}_{k,n}\) from Sect. 2.1. Define the set of population risk minimising orders by

Moreover, for any \(\psi _{k}\in \mathbb {S}_{k}\), whenever \(\textrm{E}\left\{ \left| \log f_{0}\left( X\right) \right| \right\} <\infty \), we have

and therefore minimising \(L_{k}\) over k is equivalent to minimising the Kullback–Leibler divergence \(\min _{\psi _{k}\in \mathbb {S}_{k}}\mathfrak {K}\!\left( f_{0}\Vert f_{k}\left( \cdot ;\psi _{k}\right) \right) \) over k. Equivalently,

the set of KL minimising models among \(\left\{ \mathcal{M}_{k}^{\phi }\right\} _{k\in \left[ \bar{k}\right] }\). Without any further assumptions, we have the following result.

Under A1*, A2–A4, and B1, \(\lim _{n\rightarrow \infty }\textrm{P}\!\left( \hat{k}_{n}\in \mathcal{K}^{*}\right) =1\).

Proposition 2 implies that under misspecification, any IC that satisfies B1 (and thus also any IC that satisfies B1 and B2, including the \(\epsilon \)-BIC and \(\nu \)-BIC) will eventually select an order whose best-fitting mixture model is KL-optimal among the candidate orders in \([\bar{k}]\). Besides the ICs studied in Sect. 2.4, many common ICs satisfy B1, including the AIC and the BIC. Consequently, if A1* and A2–A4 are satisfied, then all of these ICs admit the conclusion of Proposition 2.

The difference between Proposition 2 and Theorem 1 is that Proposition 2 only concludes that \(\hat{k}_{n}\) lies in the set \(\mathcal{K}^{*}\) with high probability, rather than that it converges to the smallest KL-minimising order (i.e., it is not a guarantee of parsimonious order selection). If one wishes to recover parsimonious order selection in the misspecified case, such a guarantee can be obtained via the PanIC framework of Nguyen (2024) by imposing the stronger penalty separation assumption (N-B2) from the Appendix; for instance, one may use a penalty of the form

with increasing function \(\alpha \).

We are not aware of general results proving parsimonious selection under misspecification for BIC-scale penalties, including the classical BIC, nor for our larger penalties satisfying B2, in finite mixture settings. In particular, Keribin (2000); Gassiat and Van Handel (2012), and Gassiat (2018, Sec. 4.3) all assume correct specification in their derivations. In Patilea (2001), the MLE theory of van de Geer (2000, Ch. 7) is developed for misspecified models. These results suggest a possible route to extending parsimonious order selection to misspecified settings, but carrying out such a program in the present finite-mixture framework would be nontrivial and goes beyond the scope of this work.

In this section, we provide example applications of Theorem 1 and Corollary 1 via illustrations of how the sufficient conditions are verified.

We first make our theory concrete by considering the ubiquitous Gaussian mixture model, defined by taking

with \(x\in \mathbb {X}=\mathbb {R}^{p}\) and \(\theta =\left( \mu ,\varSigma \right) \in \mathbb {T}\), where \(\mathbb {T}\) is a compact subset of \(\mathbb {R}^{p}\times \mathbb {S}_{p}^{++}\). Here, \(\mathbb {S}_{p}^{++}\) denotes the set of positive definite matrices in \(\mathbb {R}^{p\times p}\) and \(\left( \cdot \right) ^{\top }\) denotes matrix transposition. Such a compact \(\mathbb {T}\) can be constructed, as per Ritter (2014, Sec. B.23), by taking

where \(\lambda _{1}\left( \varSigma \right) \) and \(\lambda _{p}\left( \varSigma \right) \) are the smallest and largest eigenvalues of \(\varSigma \), respectively, and \(b\ge 0\) and \(c\ge 1\) are fixed constants. By this construction, A2 and A3 are naturally verified.

Next, we can write

and thus, with \(\textrm{E}\left\| X\right\| ^{2}<\infty \), we can verify A4 by taking

To verify A5*, we use standard expressions for the gradient of the normal density with respect to \(\theta \) (see, e.g., Loos 2016, Sec. 3.4.2). Write \(\theta =(\mu ,\varSigma )\), and let \(b\ge 0\) and \(c\ge 1\) be as in the definition of \(\mathbb {T}\). Define \(r(x)=\left[ \left\| x\right\| -b\right] _{+}\). Since \(\lambda _{1}(\varSigma )\ge c^{-1}\) and \(\lambda _{p}(\varSigma )\le c\), we have \(\textrm{det}(\varSigma )\ge c^{-p}\) and \(\varSigma ^{-1}\succeq c^{-1}I\), and therefore, for every \(x\in \mathbb {R}^{p}\) and \(\theta \in \mathbb {T}\),

Moreover, for the derivatives with respect to \(\mu \) and \(\varSigma \), we have

so that \(\left\| \partial \phi (x;\theta )/\partial \mu \right\| \le c\left\| x-\mu \right\| \phi (x;\theta )\). Similarly, the partial derivatives with respect to the coordinates of \(\varSigma \) can be written as \(\phi (x;\theta )\) multiplied by linear combinations of entries of \(\varSigma ^{-1}\) and of \(\varSigma ^{-1}(x-\mu )(x-\mu )^{\top }\varSigma ^{-1}\), and hence there exists a constant \(C_{0}>0\) depending only on b, c, and p, such that, for all x and \(\theta \in \mathbb {T}\),

Combining the preceding results and using \(\left\| x-\mu \right\| \le \left\| x\right\| +b\), we obtain

for a constant \(C_{1}>0\) depending only on b, c, and p. Therefore, A5* is satisfied by taking

noting that \(G_{2}\in \mathcal{L}_{1}\left( \mathfrak {m}\right) \) since it has a Gaussian tail. Since \(\textrm{E}\left\| X\right\| ^{2}<\infty \) holds under A1 for Gaussian mixture models, the conclusion of Theorem 1 holds for Gaussian mixture models whenever A1 holds. This is the archetypal example of Keribin (2000) and verifies the assumptions in the text.

The Laplace mixture model is defined by taking the density of the Laplace distribution:

with \(x\in \mathbb {R}\) and \(\theta =\left( \mu ,\gamma \right) \in \mathbb {T}\), where we can choose

which is convex and compact for each \(b>0\) and \(c>1\). Such models have been considered, for example, by Cord et al. (2006); Mitianoudis and Stathaki (2007), and Rabbani and Vafadust (2008). Here, A2 and A3 hold by construction, although it is noteworthy that (5) is a nondifferentiable function of \(\mu \).

Similarly to Sect. 3.1, we can take

to verify A4, under the condition that \(\textrm{E}\left| X\right| <\infty \). Since \(\phi \left( x;\theta \right) \) is not differentiable in \(\mu \), we seek to verify A5 instead of A5*. Write \(\theta =(\mu ,\gamma )\) and let \(b>0\) and \(c>1\) be as in the definition of \(\mathbb {T}\). Define \(r(x)=\left[ |x|-b\right] _{+}\). Then, for every \(x\in \mathbb {R}\) and \(\theta \in \mathbb {T}\), we have \(|x-\mu |\ge r(x)\) and \(\gamma \ge c^{-1}\), so that

To bound differences in \(\mu \), fix \(\gamma \in [c^{-1},c]\) and note that \(\bigl ||x-\mu _{1}|-|x-\mu _{2}|\bigr |\le |\mu _{1}-\mu _{2}|\). Using the mean value theorem for \(t\mapsto \exp \{-\gamma t\}\), we obtain

and therefore

To bound differences in \(\gamma \), fix \(\mu \in [-b,b]\) and use the mean value theorem in \(\gamma \). Since

we have

Combining the preceding bounds yields, for all \(\theta _{1},\theta _{2}\in \mathbb {T}\),

where we can take

for a suitable constant \(C_{0}>0\) depending only on b and c. Moreover, the bound on \(\phi (x;\theta )\) above shows that \(\max _{\theta \in \mathbb {T}}\phi (x;\theta )\le L_{\phi }(x)\) for \(C_{0}\) sufficiently large, and thus A5 is satisfied by taking \(G_{2}(x)=2L_{\phi }(x)\). Since \(G_{2}\in \mathcal{L}_{1}\left( \mathfrak {m}\right) \) by the exponential tails, we have the conclusion of Theorem 1 for Laplace mixture models whenever A1 is satisfied.

In this example, we observe that since \(\phi \) is not differentiable for all \(\theta \in \mathbb {T}\), the conditions (P1-b), (P2), and (P3) of Keribin (2000), defined via the partial derivatives of \(\phi \), cannot be verified.

The Student t mixture model is defined by taking the density of the Student t distribution:

with \(x\in \mathbb {R}\) and \(\theta =(\mu ,\gamma ,v)\in \mathbb {T}\), where

We can choose

which is compact for each \(b>0\), \(c>1\), and \(0<\underline{v}<\overline{v}<\infty \). Mixtures of such distributions are considered in the works of Peel and McLachlan (2000b) and Forbes and Wraith (2014), among others.

Here, A2 and A3 hold by construction, and (6) is continuously differentiable in \(\theta \) on \(\mathbb {T}\). Similarly to Sect. 3.1, we can write

and thus we can take

to verify A4, under the condition that \(\textrm{E}X^{2}<\infty \) (using \(\log (1+u)\le u\) and \(|x-\mu |\le |x|+b\)).

Since \(\phi (x;\theta )\) is differentiable in \(\theta \), we seek to verify A5\(^{*}\). Define \(r(x)=\left[ |x|-b\right] _{+}\), so that \(|x-\mu |\ge r(x)\) for all \(|\mu |\le b\). First, since \(\gamma \le c\) and \(v\in [\underline{v},\overline{v}]\), we have \(\max _{\theta \in \mathbb {T}}(c_{v}\gamma )<\infty \). Moreover, since \(\gamma \ge c^{-1}\) and \(v\le \overline{v}\), we have

Since \(v\ge \underline{v}\), the map \(t\mapsto t^{-(v+1)/2}\) is bounded above by \(t^{-(\underline{v}+1)/2}\) for \(t\ge 1\), and hence there exists a constant \(C_{0}>0\) depending only on b, c, \(\underline{v}\), and \(\overline{v}\), such that

Next, for the derivatives, we recall that for \(\theta =(\mu ,\gamma ,v)\),

and

Since \(t\mapsto t/(v+t^{2})\) is maximised at \(t=\sqrt{v}\) and \(\gamma \in [c^{-1},c]\), we have \(\gamma ^{2}|x-\mu |/\{v+\gamma ^{2}(x-\mu )^{2}\}\le c/(2\sqrt{\underline{v}})\), and thus

for constants \(C_{1},C_{2}>0\) depending only on c, \(\underline{v}\), and \(\overline{v}\). Furthermore, \(\max _{v\in [\underline{v},\overline{v}]} |\partial _{v}\log c_{v}|<\infty \) and \(\max _{v\in [\underline{v},\overline{v}]}(v+1)/(2v)<\infty \), and therefore there exists a constant \(C_{3}>0\) such that

Combining these bounds with the envelope for \(\max _{\theta }\phi (x;\theta )\) above, we obtain that

where we can take

for a suitable constant \(C_{4}>0\). Since \(\underline{v}>0\), the right-hand side has a polynomial tail of order \(|x|^{-(\underline{v}+1)}\) up to a logarithmic factor, and hence \(G_{2}\in \mathcal{L}_{1}\left( \mathfrak {m}\right) \). Therefore, A5\(^{*}\) is satisfied.

Since \(\textrm{E}X^{2}<\infty \) holds under A1 whenever the smallest degrees of freedom among the true components satisfies \(v_{0}>2\), we have the conclusion of Theorem 1 for Student t mixture models whenever A1 is satisfied with \(v_{0}>2\).

In this example, we observe that although \(\phi \) is smooth, allowing v to vary can prevent verification of the integrability requirement in (P1-b) of Keribin (2000) (see the Appendix). Let \(v_{0}\) denote the smallest degrees of freedom among the true components, so that \(f_{0}(x)\asymp |x|^{-(v_{0}+1)}\) as \(|x|\rightarrow \infty \), and fix any candidate component with degrees of freedom \(v\in [\underline{v},\overline{v}]\), for which \(\phi (x;\theta )\asymp |x|^{-(v+1)}\). Then

which diverges whenever \(v\le \frac{2}{3}v_{0}\). Moreover, for fixed v, the derivative of \(\phi (x;\theta )\) with respect to \(\gamma \) is asymptotically a nonzero constant multiple of \(\phi (x;\theta )\) as \(|x|\rightarrow \infty \), so the corresponding derivative ratio has the same integrability obstruction. Consequently, if the parameter set \(\mathbb {T}\) contains degrees of freedom \(v\le 2v_{0}/3\) (in particular, if \(\underline{v}\le 2v_{0}/3\)), then (P1-b) cannot be verified and the conditions of Keribin (2000) are not applicable in such t mixture settings. For instance, if \(v_{0}=6\) then any model class allowing \(v\le 4\) violates (P1-b).

One can extend upon finite mixture models via the mixture of regression construction which was introduced in Quandt (1972) and subsequently studied, for example, in Jones and McLachlan (1992); Wedel and DeSarbo (1995); Naik et al. (2007); Hafidi and Mkhadri (2010); Song et al. (2014); Depraetere and Vandebroek (2014), and Hui et al. (2015), among many others. A good recent account of mixture regression models appears in Yao and Xiang (2024). Suppose that we can write \(X=\left( U,Y\right) \), with \(\left( U,Y\right) :\varOmega \rightarrow \mathbb {U}\times \mathbb {Y}=\mathbb {X}\) and that we can write the PDF of X in the form

for each \(x=\left( u,y\right) \in \mathbb {X}\), where \(\rho :\mathbb {U}\times \mathbb {Y}\times \mathbb {T}\rightarrow \mathbb {R}_{\ge 0}\) characterises a conditional PDF in the sense that for every \(\left( u,\theta \right) \in \mathbb {U}\times \mathbb {T}\), \(\int _{\mathbb {Y}}\rho \left( y|u;\theta \right) \mathfrak {m}_{2}\left( \textrm{d}y\right) =1\), and \(\varphi \) is the generative PDF of U with respect to the dominating measure \(\mathfrak {m}_{1}\), where \(\mathfrak {m}=\mathfrak {m}_{1}\times \mathfrak {m}_{2}\). Note here that \(\varphi \) is not parameterised by \(\theta \), nor is it given any structure and is merely assumed to exist.

Naturally, when fitting mixture of regression models, one considers only the conditional part of (7). As such, with

for each \(k\in \left[ \bar{k}\right] \), one seeks to identify the smallest \(k_{0}\in \left[ \bar{k}\right] \), for which \(f_{0}=\mathfrak {f}_{0}\varphi \), where \(\mathfrak {f}_{0}\in \mathcal{N}_{k_{0}}^{\phi }\), or more concisely and in conforming to Sect. 1, \(f_{0}\in \mathcal{M}_{k_{0}}^{\phi }\). To this end, the usual risk, the average negative log-densities \(\ell _{k}\left( \psi _{k}\right) =-\textrm{E}\left\{ \log f_{k}\left( X;\psi _{k}\right) \right\} \), for each \(k\in \left[ \bar{k}\right] \) with \(f_{k}\in \mathcal{M}_{k}^{\phi }\), is not estimated by the natural estimator \(\ell _{k,n}\left( \psi _{k}\right) \), but instead by the empirical average negative log-conditional likelihood

where \(\mathfrak {f}_{k}\in \mathcal{N}_{k}^{\phi }\).

A useful fact is that, since \(f_{k}\left( u,y;\psi _{k}\right) =\mathfrak {f}_{k}\left( y|u;\psi _{k}\right) \varphi \left( u\right) \), we have, for each \(k\in \left[ \bar{k}\right] \) and \(\psi _{k}\in \mathbb {S}_{k}\),

Taking expectations yields \(\textrm{E}\left\{ \mathfrak {L}_{k,n}\left( \psi _{k}\right) \right\} =\ell _{k}\left( \psi _{k}\right) +\textrm{E}\left\{ \log \varphi \left( U\right) \right\} \), where the \(\textrm{E}\left\{ \log \varphi \left( U\right) \right\} \) term is the same across all k and has no influence on comparisons between risks for different values of k or different models within \(\mathcal{M}_{k}^{\phi }\) or \(\mathcal{N}_{k}^{\phi }\), since \(\varphi \) has no parameters. Moreover, the sample-level identity above shows that \(\ell _{k,n}\left( \psi _{k}\right) \) and \(\mathfrak {L}_{k,n}\left( \psi _{k}\right) \) differ by an additive term that does not depend on k or \(\psi _{k}\). Consequently, Lemma 3 and Lemma 4 remain valid with \(\mathfrak {L}_{k,n}\left( \psi _{k}\right) \) in place of \(\ell _{k,n}\left( \psi _{k}\right) \).

Furthermore, since \(P_{n}\log f=P_{n}\log \mathfrak {f}+P_{n}\log \varphi \), maximisers of the joint likelihood over \(\mathcal{M}_{k}^{\phi }\) correspond to maximisers of the conditional likelihood over \(\mathcal{N}_{k}^{\phi }\). In particular, if

then there exists

such that \(\hat{f}_{k,n}=\hat{\mathfrak {f}}_{k,n}\varphi \). If, in addition, \(\varphi \left( u\right) >0\) for all \(u\in \mathbb {U}\), then \(\log (\hat{f}_{k,n}/f_{0})=\log (\hat{\mathfrak {f}}_{k,n}/\mathfrak {f}_{0})\), and Proposition 1 is directly applicable to the maximum conditional likelihood estimator. Together, these facts imply that we can consistently estimate \(k_{0}\) by the conditional version of \(\hat{k}_{n}\):

To be precise, if \(\phi \) has form (7), and Assumptions A1–A5 hold for the corresponding joint mixture class \(\mathcal{M}_{k}^{\phi }\), then Theorem 1 holds with \(\hat{k}_{n}\) replaced by \(\tilde{k}_{n}\), with no additional assumptions beyond A1–A5. This example is general and includes many models considered in the literature, as noted at the start. In particular, under compactness restrictions on \(\mathbb {T}\) and mild integrability conditions with respect to \(\mathfrak {m}_{1}\), the verification of A4 and A5 (or A5\(^{*}\)) for regression kernels \(\rho \) proceeds exactly as in Sects. 3.2 and 3.3, with x replaced by y and with the regression location depending on u. It is possible to verify these assumptions in settings with well-behaved \(\phi \), using the results of Keribin (2000), such as in the case of normal mixture regression models (e.g., Quandt (1972) and Jones and McLachlan (1992)). However, like in the Laplace mixture example, the Laplace mixture of regression models, such as that of Song et al. (2014), have non-differentiable PDFs and thus violate the same conditions of Keribin (2000) as in the previous example. Mixture of Student t regression models such as those of Galimberti and Soffritti (2014) and Yao et al. (2014) will also incur the same problems as described in Sect. 3.3 if the degree of freedom parameters considered are too small.

Since \(\varphi \) is not parameterised by \(\theta \), it does not influence the maximisation of the likelihood over \(\psi _{k}\), and it cancels from likelihood ratios \(\log (\hat{f}_{k,n}/f_{0})\) whenever \(\varphi \left( u\right) >0\) for all \(u\in \mathbb {U}\). However, \(\varphi \) may enter the verification of A1–A5 through integrability requirements with respect to the dominating measure \(\mathfrak {m}=\mathfrak {m}_{1}\times \mathfrak {m}_{2}\). A convenient convention is to take \(\mathfrak {m}_{1}\) to be the marginal law of U, in which case \(\varphi =1\) and \(\phi \left( x;\theta \right) =\rho \left( y|u;\theta \right) \) is a density with respect to \(\mathfrak {m}=\mathfrak {m}_{1}\times \mathfrak {m}_{2}\). Under this convention, the verification of A1–A5 reduces to verifying them for \(\rho \).

Here, we provide some numerical evidence towards the comparative performance of the \(\epsilon \)-BIC and \(\nu \)-BIC versus alternatives such as the AIC, BIC, and PanIC. We consider four scenarios corresponding to our technical examples in Sects. 3.2 and 3.3, which present pathologies that do not verify the hypotheses of Keribin (2000). In each case, we simulate IID data \(\textbf{X}_{n}\) from the measure with density \(f_{0}\), a mixture with true number of components \(k_{0}\) of base density \(\phi \) with true parameter \(\psi _{k_{0}}^{*}\). Then, we estimate \(\hat{k}_{n}\) by considering mixtures of up to \(K>k_{0}\) components of the same form, using data of sizes \(n\in \left\{ 100,1000,10000\right\} \). The four scenarios considered are described in Table 2.

We repeat each scenario 200 times, and evaluate the proportions of times \(\hat{k}_{n}=k\), for each \(k\in \left[ K\right] \) and for each IC. Each IC is considered to have penalties in the form \(\textrm{pen}_{k,n}=\alpha \left( k\right) r_{n}\), where \(r_{n}\) is the shape of the penalty. For the \(\epsilon \)-BIC, we consider \(\epsilon =0.1\) and \(\epsilon =0.02\). For the \(\nu \)-BIC, we consider \(\nu =1\) (also corresponding to the \(\epsilon \)-BIC with \(\epsilon =1\)), and \(\nu =3\). For PanIC, we use the Sin–White information criterion form considered in Nguyen (2024). The shapes for the considered penalties are thus: AIC: \(r_{n}=2n^{-1}\); BIC: \(r_{n}=n^{-1}\log n\); \(\epsilon \)-BIC: \(r_{n}=n^{-1}\log ^{1+\epsilon }n\); \(\nu \)-BIC: \(r_{n}=n^{-1}\textrm{Ln}^{\circ \nu }\left( n\right) \log n\); PanIC: \(r_{n}=n^{-1/2}\sqrt{\textrm{Ln}^{\circ 2}\left( n\right) }\). The constant multiples \(\alpha \left( k\right) \) that we consider are of forms (I) \(\textrm{dim}\left( \mathbb {S}_{k}\right) /2\) (the classic BIC penalty), (II) \(\left\{ \textrm{dim}\left( \mathbb {S}_{k}\right) \right\} ^{2}/2\), (III) k (only the number of mixture components), and (IV) \(k^{2}\). The results are recorded in Tables 3, 4, 5, 6.

From Tables 3–6, we observe, as expected, that as n increases, the proportions of times that the \(\epsilon \)-BIC, \(\nu \)-BIC, and PanIC estimate the correct number of components increase, and typically become perfectly accurate, except in cases where the penalty and shape combinations are too severe and thus consistency is not realised in practice even at \(n=10000\), such as for PanIC with Type II penalties or the \(\nu =1\) case of the \(\nu \)-BIC with Type II penalties. We observe in our experiments that the Type III penalties when applied with the \(\epsilon \)-BIC or \(\nu \)-BIC shapes tend to be too small to realise perfect accuracy at \(n=10000\), although it is clear that the proportion of times \(\hat{k}_{n}=k_{0}\) is increasing in each case. In terms of good choices for \(\alpha \), we observe that Type I and Type IV forms appear to work well with all penalty shapes, with Type II often being too large and Type III being too small.

Perhaps the most surprising outcome is that for some configurations, the AIC and the BIC can produce perfect or near perfect estimation of \(k_{0}\) for large n. Indeed, in the case of the AIC, as we shall show in Proposition 7 (albeit in a simplified setting that does not include the current experimental scenarios), it is impossible to guarantee that the proportion of accurate estimates of \(k_{0}\) goes to one. This points to either this being a small sample phenomenon and thus the true proportion of overestimating \(k_{0}\) is possibly small but non-zero, or that there exist settings where the AIC is consistent, although we believe that the former is more likely than the latter given our results in Sect. 4.2. In the case of the BIC, we argue that it is unknown whether the conditions of Keribin (2000) are necessary and sufficient, and that there is potentially a gap between our theory and the theory of Keribin (2000). Furthermore, if we observe the results of the \(\epsilon \)-BIC with \(\epsilon =0.02\) and the \(\nu \)-BIC with \(\nu =3\), we note that the outcomes are either exactly the same as or very close to those of the BIC for all n, due to the fact that, numerically, the criteria are only infinitesimally different for any practical size n. We, however, conclude that there is indeed evidence that the infinitesimal variations from the BIC, via the \(\epsilon \)-BIC and \(\nu \)-BIC constructions, are a harmless and practically costless way to guarantee consistency in settings where such a guarantee may not otherwise be afforded by the results of Keribin (2000). Overall, these simulations lead to our recommendation, at least in these simulation scenarios, that the \(\epsilon \)-BIC with \(\epsilon =0.02\) and the \(\nu \)-BIC with \(\nu =3\), together with \(\alpha \) of Type I or Type IV, provide good finite sample performance while enabling the guarantees of Corollary 1.

Theorem 1 provides a very general and flexible approach to order selection for mixture models, and as pointed out in Remarks 4 and 5, there are many possible choices for penalty calibration beyond our proposed \(\nu \)-BIC and \(\epsilon \)-BIC proposals. A natural question to ask is whether there are optimal choices for the penalty. In particular, from a minimax perspective (cf. Gassiat and Van Handel (2012) and Gassiat (2018), Ch. 4), we would like to determine if there is a smallest penalty that satisfies B1 and B2. Such a penalty would provide model selection consistency while minimally modifying the maximum likelihood as the axis for model choice.

Unfortunately, while such a minimal penalty is desirable, we can prove that it does not exist. In particular, for fixed \(\bar{k}\ge 2\), let \(\mathscr {P}\) denote the collection of penalty arrays \(\textrm{pen}=\left\{ \textrm{pen}_{k,n}:k\in \left[ \bar{k}\right] ,n\in \mathbb {N}\right\} \) that satisfy B1 and B2. We endow \(\mathscr {P}\) with the domination order: \(\textrm{pen}^{\left( 1\right) }\preceq \textrm{pen}^{\left( 2\right) }\) if and only if there exists an \(n_{0}\in \mathbb {N}\) such that for all \(n\ge n_{0}\) and for all \(k\in \left[ \bar{k}\right] \), \(\textrm{pen}_{k,n}^{\left( 1\right) }\le \textrm{pen}_{k,n}^{\left( 2\right) }\). Call \(\textrm{pen}^{*}\in \mathscr {P}\) universally optimal if it is a least element of \(\left( \mathscr {P},\preceq \right) \); i.e., \(\textrm{pen}^{*}\preceq \textrm{pen}\), for each \(\textrm{pen}\in \mathscr {P}\).

The partial order \(\left( \mathscr {P},\preceq \right) \) has no minimal element.

More precisely, for every \(\textrm{pen}\in \mathscr {P}\), there exists a \(\widetilde{\textrm{pen}}\in \mathscr {P}\) such that \(\widetilde{\textrm{pen}}\preceq \textrm{pen}\) and such that the domination is strict for at least one order for all sufficiently large n; indeed, it is strict for every \(k\ge 2\) for all sufficiently large n. In particular, this shows that there is no optimal choice of penalty rate and \(\alpha \) function among the class of IC penalties characterised by B1 and B2 alone, and any penalty that satisfies B1 and B2 admits an asymptotically smaller penalty that also satisfies B1 and B2.

We may hope that if we have some fixed \(\alpha :\mathbb {N}\rightarrow \mathbb {R}_{>0}\) strictly increasing, and we define our penalty in the form \(\textrm{pen}_{k,n}^{\alpha }=\alpha \left( k\right) r_{n}\), where \(\left( r_{n}\right) _{n}\subset \mathbb {R}_{>0}\), then there may be some minimal choice of \(\left( r_{n}\right) _{n}\). Unfortunately, this also is not true as per the following result.

For penalties of the form \(\textrm{pen}_{k,n}^{\alpha }=\alpha \left( k\right) r_{n}\), if B1 and B2 hold then necessarily \(r_{n}\rightarrow 0\) and \(\left( n/\log n\right) r_{n}\rightarrow \infty \), as \(n\rightarrow \infty \). Moreover, for each such \(\left( r_{n}\right) _{n}\), there exists another sequence \(\left( \tilde{r}_{n}\right) _{n}\) such that \(\tilde{r}_{n}=o\left( r_{n}\right) \), while still satisfying \(\tilde{r}_{n}\rightarrow 0\) and \(\left( n/\log n\right) \tilde{r}_{n}\rightarrow \infty \), as \(n\rightarrow \infty \).

Let \(a_{n}=\left( n/\log n\right) r_{n}\rightarrow \infty \) and set \(\tilde{r}_{n}=r_{n}/\sqrt{a_{n}}\). Then \(\tilde{r}_{n}/r_{n}=a_{n}^{-1/2}\rightarrow 0\), hence \(\tilde{r}_{n}=o\left( r_{n}\right) \) and \(\tilde{r}_{n}\rightarrow 0\). Lastly, \(\left( n/\log n\right) \tilde{r}_{n}=\sqrt{a_{n}}\rightarrow \infty \), satisfying B2, since \(\alpha \left( l\right) -\alpha \left( k\right) >0\) for each \(l>k\), as required. \(\square \)

Indeed, the same penalty flexibility issue is not only a feature of our consistency theorems but also of the PanIC approach of Nguyen (2024) and Westerhout et al. (2024), the preceding results of Sin and White (1996), and even the mixture model BIC consistency results of Keribin (2000); Gassiat and Van Handel (2012), and Gassiat (2018). This is also a problem in the finite sample model selection approaches of Birgé and Massart (2007) and Massart (2007) (implemented in the mixture model literature in the works of Maugis and Michel (2011); Maugis-Rabusseau and Michel (2013), and Nguyen et al. (2022)). In the finite sample setting, some methods are available for empirical calibration of penalties, such as via the so-called slope heuristic (cf. Baudry et al. (2012) and Arlot (2019)). Although there is some evidence that these methods work well in empirical studies, they provide no technical or concrete guarantees in the mixture model setting.

In Yang (2005), it was demonstrated that, in the regression context, there is a fundamental tension between the AIC and consistent estimators such as the BIC. In particular, it was shown that in the context of choosing between nested regression models, the chosen model selected via the AIC is minimax optimal in the mean squared error sense, while it is impossible for any fitted model chosen via a model consistent estimator to achieve the same minimax optimality. In this subsection, we will prove that the same phenomenon arises in the mixture model order selection setting.

To proceed, we will consider a simplified version of the example in Sect. 3.1. In particular, let

be the density of the univariate normal law with unit variance \(\textrm{N}\left( \mu ,1\right) \). For fixed \(b>0\) and \(\underline{\pi }\in \left( 0,1/2\right) \), let \(\mathcal{F}=\mathcal{M}_{1}\cup \mathcal{M}_{2}\), where \(\mathcal{M}_{1}=\left\{ \phi \left( \cdot ;\mu \right) :\mu \in \left[ -b,b\right] \right\} \), and

Let \(\textbf{X}_{n}\) consist of IID replicates from the law with density (with respect to the Lebesgue measure) \(f_{0}\in \mathcal{F}\). As in Sect. 3.4, we let \(\hat{f}_{k,n}\in \arg \max _{f\in \mathcal{M}_{k}}P_{n}\log f\) for \(k\in \left\{ 1,2\right\} \), and with \(f\mapsto \ell _{n}\left( f\right) =-n^{-1}\sum _{i=1}^{n}\log f\left( X_{i}\right) \), we write

Here, we define AIC-like penalties via \(\textrm{pen}_{k,n}^{\textrm{AIC}}=a\textrm{dim}\left( \mathbb {S}_{k}\right) n^{-1}\), where \(a>0\) is sufficiently large, and \(\textrm{dim}\left( \mathbb {S}_{1}\right) =1\) and \(\textrm{dim}\left( \mathbb {S}_{2}\right) =3\), in this case. As per the usual treatment on the topic (see, e.g., Wainwright (2019), Ch. 15), we define the squared-Hellinger minimax risk as

where \(\textrm{E}_{f_{0}}\) denotes the expectation over \(\textbf{X}_{n}\) with respect to the law with density \(f_{0}\). The infimum is taken over the class of all measurable estimators \(\tilde{f}_{n}\). The following result states that the AIC-like estimator \(\hat{f}_{n}^{\textrm{AIC}}\) achieves the minimax risk.

For each sufficiently large \(a>0\), there exist constants \(0<c<C<\infty \) depending on a, b, and \(\underline{\pi }\) such that, for all sufficiently large n,

In contrast to the AIC-like penalties, we note that the BIC, and more generally any penalty satisfying B1 together with B2, verifies the condition

where \(\varDelta _{n}=\textrm{pen}_{2,n}-\textrm{pen}_{1,n}\). The same implication holds for the average-risk form of Condition (C1) of Keribin (2000) and for penalties satisfying (N-B1)–(N-B2) of Nguyen (2024). With \(\tilde{f}_{n}=\hat{f}_{\hat{k}_{n},n}\), where \(\hat{k}_{n}=\arg \min _{k\in \left\{ 1,2\right\} }\left\{ \ell _{n}\left( \hat{f}_{k,n}\right) +\textrm{pen}_{k,n}\right\} \), we have the following result.

If (8) holds then there exists a constant \(c>0\) depending only on b and \(\underline{\pi }\), such that for all sufficiently large n,

Proposition 6 implies that any selector \(\hat{k}_{n}\) that satisfies (8) cannot achieve the minimax Hellinger risk over \(\mathcal{F}\). Indeed, from Proposition 5, we have that \(\mathscr {R}_{n}^{*}\asymp n^{-1}\), which, combined with Proposition 6, yields

Thus, the worst-case risk of \(\tilde{f}_{n}\) is asymptotically larger than the minimax risk by an unbounded factor. A natural follow-up question is whether it is possible to construct a consistent selector that does not verify condition (8). The following result demonstrates that this is impossible in the current setting.

Suppose that \(\hat{k}_{n}\) is the penalised likelihood selector on \(\mathcal{F}=\mathcal{M}_{1}\cup \mathcal{M}_{2}\) with penalty gap \(\varDelta _{n}=\textrm{pen}_{2,n}-\textrm{pen}_{1,n}\). If \(\hat{k}_{n}\) is order-consistent in the sense that, for every \(f_{0}\in \mathcal{M}_{1}\), \(\textrm{P}_{f_{0}}(\hat{k}_{n}=2)\rightarrow 0\), and, for every \(f_{0}\in \mathcal{M}_{2}\setminus \mathcal{M}_{1}\), \(\textrm{P}_{f_{0}}(\hat{k}_{n}=1)\rightarrow 0\), then (8) necessarily holds.

The IC approach to model selection and order selection is ubiquitous across many domains of statistical modeling, with particular popularity and applicability in the context of finite mixture models. Among the IC approaches, the BIC is often the most commonly used, and in the mixture model setting where the maximum model size \(\bar{k}\) is known, its consistency has been proved in Keribin (2000). In this work, we substantially relax the assumptions made by Keribin (2000) via minor modifications of the BIC, to produce the \(\nu \)-BIC and \(\epsilon \)-BIC. These ICs minimally modify the penalty term of the BIC but allow for much broader applicability, notably being applicable without any differentiability and requiring substantially weaker regularity than high-order differentiability and moment conditions. Via example applications, we show that these ICs are consistent estimators in the popular setting of Gaussian mixture models, and also in settings that fall outside the scope of Keribin (2000), including non-differentiable Laplace mixture models, heavy-tailed t-mixture models, and mixtures of regression models. In addition, we provide a complementary result under misspecification: when \(f_{0}\) need not belong to any \(\mathcal{M}_{k}^{\phi }\), any IC satisfying B1 will eventually select an order whose best-fitting mixture model is Kullback–Leibler optimal among the candidate orders in \([\bar{k}]\) (Proposition 2).

As we have discussed in Sect. 1, for suitable and reasonable choices of \(\nu \) and \(\epsilon \), the respective \(\nu \)-BIC and \(\epsilon \)-BIC are numerically indistinguishable from the BIC in practical settings. We accompany this observation with numerical simulations in Sect. 3, which compare the \(\epsilon \)-BIC and \(\nu \)-BIC against alternatives such as the AIC, BIC, and PanIC in scenarios that exhibit the pathologies highlighted in our technical examples. In these experiments, we observe that choices such as \(\epsilon =0.02\) and \(\nu =3\) are, for all practical values of n, effectively the same as the BIC numerically, while still inheriting the theoretical guarantees of Corollary 1 in settings where the hypotheses of Keribin (2000) are not available. More generally, the simulation study provides some empirical guidance on the practical calibration of \(\alpha \) and on the behaviour of competing criteria in difficult regimes.

Given the potentially small numerical difference between the \(\nu \)-BIC and \(\epsilon \)-BIC and the BIC, we do not seek to push the use of either criterion over the BIC, which many practitioners are already accustomed to and use regularly. However, we wish to position our results as technical tools for explaining why penalties that are numerically close to the BIC can remain valid even when the stringent assumptions of Keribin (2000) are not met. For example, one may use our theory to justify BIC-like modifications for mixtures with non-differentiable components \(\phi \), such as triangular mixtures (Karlis and Xekalaki, 2008; Nguyen and McLachlan, 2016) and Laplace-based mixture models (Franczak et al., 2013; Song et al., 2014; Azam and Bouguila, 2020), as well as mixtures with components \(\phi \) with limited regularity, such as t-distribution based models (Peel and McLachlan, 2000b; Lin et al., 2007; Lo and Gottardo, 2012; Forbes and Wraith, 2014; Yao et al., 2014). In the misspecified case, Proposition 2 further indicates that, under mild conditions, common ICs (including the AIC and the BIC) still target Kullback–Leibler optimal approximating models, albeit without a guarantee of parsimonious order selection.

Indeed, via Theorem 1, it is possible to construct multitudes of ICs that exhibit consistency in mixture model settings, but the finite sample performance of any such penalty cannot be deduced from our general consistency theory alone. In particular, as formalised in Sect. 4.2, there are two structural limitations that consistency results do not resolve. First, under the natural domination order, there is no universally minimal penalty array satisfying B1 and B2, so consistency does not single out an optimal BIC-scale penalty calibration. Second, we establish a tension between order consistency and minimax optimality in mixtures like that of Yang (2005): in a simple Gaussian mixture family, an AIC-like criterion achieves the parametric minimax Hellinger risk, whereas any order-consistent criterion must be minimax-inefficient by an unbounded factor, typically logarithmic for BIC-scale penalties. To the best of our knowledge, this AIC/BIC tension has not previously been stated explicitly in the finite-mixture order selection setting. In the finite sample setting, some methods are available for empirical calibration of penalties, such as via the so-called slope heuristic (cf. Baudry et al. (2012) and Arlot (2019)), although such methods provide no general guarantees in the mixture model setting.

Although the assumptions of Theorem 1 are weaker than those of Keribin (2000), it is still possible to deduce consistent model selection under further relaxations of conditions. For example, even without continuity of the component density \(\phi \) or independence of the data \(\textbf{X}_{n}\), or when the loss function is not necessarily the negative log-density, we can deduce consistency using the theory of Westerhout et al. (2024), at the expense of increasing the rate of the penalisation term to \(\tilde{O}\left( n^{-1/2}\right) \). An immediate direction of study is to understand whether further relaxations of the assumptions are possible without such a pronounced increase in the penalisation rate. Another natural direction, motivated by truth-changing and locally misspecified frameworks in modern model selection (for example the focused information criterion literature of Claeskens and Hjort (2003); Lohmeyer et al. (2019); Pandhare and Ramanathan (2020a, 2020b) and related high-dimensional settings such as Gao and Carroll (2017)), is to extend our results to triangular-array regimes where the effective data-generating mechanism, or the KL-optimal order, may vary with n. We expect such extensions to be possible, but they would require triangular-array analogues of the empirical process arguments underpinning our uniform convergence and overfitting rate results. In addition, under misspecification, Proposition 2 only ensures that \(\hat{k}_{n}\) eventually lies in the set of KL-minimising orders, and we are not aware of general results proving parsimonious order selection under misspecification for BIC-scale penalties. It would therefore be of interest to understand whether parsimonious selection is achievable under misspecification without moving to larger penalties, or whether such a limitation is inherent.

Another immediate extension to our current work is to obtain the corresponding regularity conditions for consistent order selection of mixture of experts models (cf. Jacobs et al. (1991); Yuksel et al. (2012), and Nguyen and Chamroukhi (2018)), which extend upon the finite mixture models by also parameterising the mixing coefficients \(\pi _{1},\dots ,\pi _{k}\). Such an extension would add to the growing body of theoretical results regarding model selection and model uncertainty in mixture of experts models, including the results of Khalili (2010); Nguyen et al. (2022, 2023c, 2024c, 2023, 2024a, 2024b); Thai et al. (2025), and Westerhout et al. (2024), for example. Lastly, a limitation of our current results is the requirement that some maximum order \(\bar{k}\) be specified. Indeed such an assumption can be omitted at the great expense of strong control of local bracketing entropy of normalisations of the mixture model classes, as in Gassiat and Van Handel (2012) and Gassiat (2018, Sec. 4.3). It remains open as to whether weaker assumptions are available, even for penalties with the larger rate \(\tilde{O}\left( n^{-1/2}\right) \).