Testing the Order of Multivariate Normal Mixture Models

The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.

105,739 characters

Testing the Order of Multivariate Normal Mixture Models Hiroyuki Kasahara Vancouver School of Economics University of British Columbia [email removed] Katsumi Shimotsu Faculty of Economics University of Tokyo [email removed]


\bibliographystyle{asa}

{
\singlespacing
\title{\textbf{Testing the Order of Multivariate Normal Mixture Models}
\author{Hiroyuki Kasahara\thanks{Address for correspondence: Hiroyuki Kasahara, Vancouver School of Economics, University of British Columbia, 997-1873 East Mall, Vancouver, BC V6T 1Z1, Canada.   The authors thank Pengfei Li and participants at the conference on Advances in Finite Mixture and other Non-regular Models, Guilin, China in 2018  for helpful comments and the Institute of Statistical Mathematics for the facilities and the use of SGI ICE X. This research is support by the Natural Science and Engineering Research Council of Canada and JSPS Grant-in-Aid for Scientific Research (C) No. 26380267.  }\\
Vancouver School of Economics\\
University of British Columbia\\
[email removed] \and Katsumi Shimotsu\\
Faculty of Economics \\
University of Tokyo\\
[email removed]
}}
\maketitle
}
\begin{abstract}
Finite mixtures of multivariate normal distributions have been widely used in empirical applications in diverse fields such as statistical genetics and statistical finance. Testing the number of components in multivariate normal mixture models is a long-standing challenge even in the most important case of testing homogeneity. This paper develops likelihood-based tests of the null hypothesis of $M_0$ components against the alternative hypothesis of $M_0 + 1$ components for a general $M_0 \geq 1$. For heteroscedastic normal mixtures, we propose an EM test and derive the asymptotic distribution of the EM test statistic. For homoscedastic normal mixtures, we derive the asymptotic distribution of the likelihood ratio test statistic. We also derive  the asymptotic distribution of the likelihood ratio test statistic and EM test statistic under local alternatives and show the validity of parametric bootstrap. The simulations show that the proposed test has good finite sample size and power properties.
\end{abstract}

Key words: asymptotic distribution; EM test; likelihood ratio test; multivariate normal mixture models; number of components

\section{Introduction}


Finite mixtures of multivariate normal distributions have been widely used in empirical applications in diverse fields such as statistical genetics and statistical finance. Comprehensive surveys on theoretical properties and applications can be found, for example, \citet{lindsay95book}, \citet{mclachlanpeel00book}, and \citet{fruehwirth06book}.

In many applications of finite mixture models, the number of components is of substantial interest. In multivariate normal mixture models, however, testing for the number of components has been an unsolved problem even in the most important case of testing homogeneity. For general finite mixture models, the asymptotic distribution of the likelihood ratio test statistic (LRTS) has been derived as a functional of the Gaussian process \citep{dacunha99as, liushao03as, zhuzhang04jrssb, azais09esaim}. These results are not applicable to normal mixtures because normal mixtures have an undesirable mathematical property that invalidates key assumptions in these works \citep{chenlifu12jasa}. In particular, the normal density with mean $\mu$ and variance $\sigma^2$, $f(y;\mu,\sigma^2)$, has the property $\frac{\partial^2}{\partial \mu \partial \mu} f(y;\mu,\sigma^2) = 2\frac{\partial}{\partial \sigma^2} f(y;\mu,\sigma^2)$. This leads to the loss of ``strong identifiability'' condition introduced by \citet{chen95as}. As a result, neither Assumption (P1) of \citet{dacunha99as} nor Assumption 7 of \citet{azais09esaim} holds, and Assumption 3 of \citet{zhuzhang04jrssb} is violated, while Corollary 4.1 of \citet{liushao03as} does not hold in normal mixtures. Heteroscedastic normal mixture models have an additional problem called the infinite Fisher information problem \citep{lcm09bm} that the score of the LRTS has infinite variance if the range of the variance is unrestricted.

This paper develops likelihood-based tests of the null hypothesis of $M_0$ components against the alternative hypothesis of $M_{0} + 1$ components for a general $M_0 \geq 1$ in multivariate normal mixtures. We consider both heteroscedastic and homoscedastic mixtures. For heteroscedastic normal mixtures, we propose an EM test by building on the EM approach pioneered by \citet{lcm09bm} and \citet{lichen10jasa}. The asymptotic null distribution of the proposed EM test statistic is shown to be the maximum of $M_0$ random variables, each of which is a projection of a Gaussian random variable on a cone. For homoscedastic normal mixtures, we derive the asymptotic distribution of the LRTS because homoscedastic normal mixtures do not suffer from the infinite Fisher information problem.

In univariate heteroscedastic normal mixtures, \citet{chenli09as} develop an EM test for $M_0 = 1$ against $M_0=2$, and \citet{chenlifu12jasa} develop an EM test for testing $H_0:M = M_0$ against $H_A:M > M_0$. Our result may be viewed as generalization of \citet{chenli09as} to the multivariate case. \citet{kasaharashimotsu15jasa} develop an EM test for testing $H_0:M = M_0$ against $H_A:M = M_0+1$ for general $M_0 \geq 1$ in finite normal mixture regression models. In univariate homoscedastic normal mixtures, \citet{chenchen03sinica} derive the asymptotic distribution of the LRTS. Our results generalize the results in \citet{chenchen03sinica} to multivariate homoscedastic normal mixtures. For some specific models such as binomial mixtures, the asymptotic distribution of the LRTS has been derived by, for example,  \citet{ghoshsen85book, chernofflander95jspi, lemdanipons97spl, chenchen01cjstat, chenchen03sinica, cck04jrssb, garel01jspi, garel05jspi}.

The remainder of this paper is organized as follows. Section 2 introduces the likelihood ratio test for heteroscedastic multivariate normal mixture models as a precursor of the EM test and derives the asymptotic distribution of the LRTS. Section 3 introduces the EM test and derives the asymptotic distribution of the EM test statistic. Section 4 derives the asymptotic distribution of the LRTS and EM test statistics under local alternatives and Section 5 shows the validity of parametric bootstrap. Section 6 analyzes  homoscedastic multivariate normal mixture models. Section 7 reports the simulation results and provides empirical applications. Appendix A contain proofs, and Appendices B--D collect auxiliary results.

We collect notation. Let $:=$ denote ``equals by definition.'' Boldface letters denote vectors or matrices. For a matrix $\boldsymbol{B}$, denote its $(i,j)$ element by $B_{ij}$, and let $\lambda_{\min}(\boldsymbol{B}) $ and $\lambda_{\max}(\boldsymbol{B})$ be the smallest and the largest eigenvalue of $\boldsymbol{B}$, respectively. For a $k$-dimensional vector $\boldsymbol{x} = (x_1,\ldots,x_k)^{\top}$ and a matrix $\boldsymbol{B}$, define $|\boldsymbol{x}| := (\boldsymbol{x}^{\top}\boldsymbol{x})^{1/2}$ and $|\boldsymbol{B}| := (\lambda_{\max}(\boldsymbol{B}^{\top}\boldsymbol{B}))^{1/2}$. Let $\boldsymbol{x}^{\otimes k} := \boldsymbol{x} \otimes \boldsymbol{x} \otimes \cdots \otimes \boldsymbol{x}$ ($k$ times). Let $\mathbb{I}\{A\}$ denote an indicator function that takes value 1 when $A$ is true and 0 otherwise. $\mathcal{C}$ denotes a generic nonnegative finite constant whose value may change from one expression to another. Given a sequence $\{f(\boldsymbol{Y}_i)\}_{i=1}^n$, let $\nu_n(f(\boldsymbol{y})) := n^{-1/2} \sum_{i=1}^n [f(\boldsymbol{Y}_i) - Ef(\boldsymbol{Y}_i)]$ and $P_n(f(\boldsymbol{y})) := n^{-1} \sum_{i=1}^n f(\boldsymbol{Y}_i)$. All the limits are taken as $n \to \infty$ unless stated otherwise.

\section{Heteroscedastic multivariate  finite normal mixture models}\label{sec:hetero}

Denote the density of a $d$-variate normal distribution with mean $\boldsymbol{\mu}+ \boldsymbol{\gamma}^\top\boldsymbol{z}$ and variance $\boldsymbol{\Sigma}$ by
\begin{equation} \label{normal_density}
f(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\mu},\boldsymbol{\Sigma}):= (2\pi)^{-\frac{d}{2}}(\det \boldsymbol{\Sigma})^{-\frac{1}{2}} \exp \left(-\frac{(\boldsymbol{x}-\boldsymbol{\mu}-\boldsymbol{\gamma}^\top\boldsymbol{z})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}-\boldsymbol{\gamma}^\top\boldsymbol{z})}{2}\right),
\end{equation}
where $\boldsymbol{x}$ and $\boldsymbol{\mu}$ are $d \times 1$, $\boldsymbol{\gamma}$ is $d\times p$, and  $\boldsymbol{z}$ is $p \times 1$. Let $\Theta_{\boldsymbol{\gamma}} \subset \mathbb{R}^{dp}$, $\Theta_{{\boldsymbol{\mu}}} \subset \mathbb{R}^{d}$, and $\Theta_{\boldsymbol\Sigma} \subset \mathbb{S}^d_+$ denote the space of $\boldsymbol{\gamma}$, ${\boldsymbol{\mu}}$, and $\boldsymbol{\Sigma}$, respectively, where $\mathbb{S}^d_+$ denotes the space of $d\times d$ positive definite matrices. For $M \geq 2$, denote the density of $M$-component finite normal mixture distribution as:
\begin{equation}
f_M(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_M) : = \sum_{j = 1}^M \alpha_j f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma},{\boldsymbol{\mu}}_j,\boldsymbol{\Sigma}_j),\label{general_model}
\end{equation}
where $\boldsymbol{\vartheta}_M := ({\boldsymbol{\alpha}},\boldsymbol{\gamma},{\boldsymbol{\mu}}_1,\ldots,{\boldsymbol{\mu}}_{M},\boldsymbol{\Sigma}_1,\ldots,\boldsymbol{\Sigma}_M)$ with ${\boldsymbol{\alpha}} : = (\alpha_1,\ldots,\alpha_{M-1})^{\top}$, and $\alpha_{M}$ being determined by $\alpha_{M}: = 1-\sum_{j = 1}^{M - 1} \alpha_j$. $\boldsymbol{\mu}_j$ and $\boldsymbol{\Sigma}_j$ are mixing parameters that characterize the $j$-th component, and $\alpha_j$s are mixing probabilities. $\boldsymbol{\gamma}$ is the coefficient of the covariate $\boldsymbol{z}$, and $\boldsymbol{\gamma}$ is assumed to be common to all the components. Define the set of admissible values of ${\boldsymbol{\alpha}}$ by $\Theta_{{\boldsymbol{\alpha}}} :=\{{\boldsymbol{\alpha}}: \alpha_j \geq 0, \sum_{j = 1}^{M - 1} \alpha_j \in [0,1] \}$, and let the space of ${\boldsymbol{\vartheta}}_M$ be $\Theta_{{\boldsymbol{\vartheta}}_M} := \Theta_{{\boldsymbol{\alpha}}}\times \Theta_{\boldsymbol{\gamma}}\times\Theta_{{\boldsymbol{\mu}}}^M\times\Theta_{\Sigma}^M$.

The number of components $M$ is the smallest number such that the data density admits the representation (\ref{general_model}). Our objective is to test
\[
H_0:\ M = M_0\quad \text{against}\quad H_A: M = M_0 + 1.
\]


\subsection{Likelihood ratio test of $H_0: M = 1$ against $H_A: M = 2$} \label{section:homogeneity}

As a precursor of the EM test developed in Section \ref{section:emtest}, this section establishes the asymptotic distribution of the LRTS for testing the null hypothesis $H_0: M = 1$ against $H_A: M = 2$ when the data are from $H_0$.


We consider a random sample of $n$ independent observations $\{\boldsymbol{X}_i,\boldsymbol{Z}_i\}_{i = 1}^n$ from the true one-component density $f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma}^*, \boldsymbol{\mu}^*,\boldsymbol{\Sigma}^{*})$. Here, the superscript $*$ signifies the true parameter value. Let a two-component mixture density with ${\boldsymbol{\vartheta}}_2 = (\alpha,\boldsymbol{\gamma}, \boldsymbol{\mu}_1, \boldsymbol{\mu}_2,\boldsymbol{\Sigma}_1,\boldsymbol{\Sigma}_2) \in \Theta_{{\boldsymbol{\vartheta}}_2}$ be
\begin{equation}
f_2(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\vartheta}_2) := \alpha f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma} ,\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1) + (1-\alpha) f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma},\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2). \label{two-component}
\end{equation}
We partition the null hypothesis $H_0: m = 1$ into two as follows:
\[
H_{01}: (\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_2)\ \text{ and }\ H_{02}: \alpha(1-\alpha) = 0.
\]
In the following, we focus on testing $H_{01}:(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_2)$ because, as discussed in \citet{chenli09as}, the Fisher information for testing $H_{02}$ is not finite unless the range of $\det(\boldsymbol{\Sigma}_1)/\det(\boldsymbol{\Sigma}_2)$ is restricted.

The log-likelihood function for testing $H_{01}:(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_2)$ is unbounded if $\det(\boldsymbol{\Sigma}_1)$ and $\det(\boldsymbol{\Sigma}_1)$ are not bounded away from 0 \citep{hartigan85book}. Therefore, we consider a maximum penalized likelihood estimator (PMLE) introduced by \citet{chentan09jmva}. Similar to \citet{chentan09jmva}, we use the following penalty function with $M=2$:
\begin{equation} \label{pen_pmle}
p_n(\boldsymbol{\vartheta}_M) = \sum_{m=1}^M p_{n}(\boldsymbol{\Sigma}_m;\widehat{\boldsymbol{\Omega}}) = \sum_{m=1}^M -a_n\left\{ \text{tr}(\widehat{\boldsymbol{\Omega}}\boldsymbol{\Sigma}_m^{-1}) - \log ( \det(\widehat{\boldsymbol{\Omega}}\boldsymbol{\Sigma}_m^{-1})) -d \right\},
\end{equation}
where $\widehat{\boldsymbol{\Omega}}$ is the maximum likelihood estimator (MLE) of $\boldsymbol{\Sigma}$ from the one-component model, and $a_n$ is a non-random sequence such that $a_n \geq 1/n$ and $a_n =o(n)$.  Let $\widehat{\boldsymbol{\vartheta}}_2$ denote the PMLE that maximizes $PL_n({\boldsymbol{\vartheta}}_2):= \sum_{i = 1}^n f_2(\boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_2) + p_n(\boldsymbol{\vartheta}_2)$.

\begin{assumption} \label{assn_consis} $\boldsymbol{Z}$ has finite second moment, and $\Pr(\boldsymbol{\gamma}^\top \boldsymbol{Z}_i \neq \boldsymbol{\gamma}^{*\top} \boldsymbol{Z}_i )>0$ for any $\boldsymbol{\gamma} \neq \boldsymbol{\gamma}^*$.
\end{assumption}
Model (\ref{two-component}) yields the true density $f(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}^*,\boldsymbol{\Sigma}^{*})$ if ${\boldsymbol{\vartheta}}_2$ lies in the set $\Theta_2^* := \{{\boldsymbol{\vartheta}}_2 \in \Theta_{{\boldsymbol{\vartheta}}_2}: \{ (\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2) = (\boldsymbol{\mu}^*,\boldsymbol{\Sigma}^{*}), \boldsymbol{\gamma}=\boldsymbol{\gamma}^*\}\ \text{or}\ \{\alpha =1, (\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}^*,\boldsymbol{\Sigma}^{*}), \boldsymbol{\gamma}=\boldsymbol{\gamma}^*\}\ \text{or}\ \{\alpha=0, (\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2) = (\boldsymbol{\mu}^*,\boldsymbol{\Sigma}^{*}), \boldsymbol{\gamma}=\boldsymbol{\gamma}^*\}\}$. The following proposition shows the consistency of $\widehat{\boldsymbol{\vartheta}}_2$.
\begin{proposition} \label{P-consis}
Suppose that Assumption \ref{assn_consis} holds. Then, under the null hypothesis $H_0: M=1$, $\inf_{{\boldsymbol{\vartheta}}_2 \in \Theta_2^*} |\widehat {\boldsymbol{\vartheta}}_2 - {\boldsymbol{\vartheta}}_2| \rightarrow_p 0$.
\end{proposition}
In testing $H_{01}$, the standard asymptotic analysis of the LRTS breaks down because the Fisher information matrix is degenerate. This is due to the fact that, for any $\bar {\boldsymbol{\vartheta}}_2$ such that $(\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1)=(\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2)$, the derivatives of the density of different orders  are linearly dependent as
\begin{align*}
&\nabla_{\boldsymbol{\mu}_1} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2) =\frac{\alpha}{1-\alpha} \nabla_{\boldsymbol{\mu}_2} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2),\ \nabla_{\boldsymbol{\Sigma}_1}f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2) = \frac{\alpha}{1-\alpha} \nabla_{\boldsymbol{\Sigma}_2}l(y|\boldsymbol{x,z};\bar {\boldsymbol{\vartheta}}_2), \nonumber \\
&\nabla_{\mu_{1i}\mu_{1j}} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2) = 2\nabla_{\Sigma_{1,ij}} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2), \quad \nabla_{\mu_{2i}\mu_{2j}} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2) = 2\nabla_{\Sigma_{2,ij}} f_2(\boldsymbol{x}|\boldsymbol{z};\bar{\boldsymbol{\vartheta}}_2).
\end{align*}
This dependence leads to the loss of strong identifiability and causes substantial difficulties in existing literature.

We analyze the LRTS for testing $H_{01}: (\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_1)=(\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_2)$ by developing a higher-order approximation of the log-likelihood function  through an ingenious reparameterization that  extends the result of \citet{rotnitzky00bernoulli} and \citet{kasaharashimotsu15jasa}. Collect the unique elements in $\boldsymbol{\Sigma}$ into a $d(d+1)/2$-vector
\begin{equation*}
\begin{aligned}
\boldsymbol{v} &= (v_{11},v_{12},\ldots,v_{1d},v_{22},v_{23},\ldots,v_{2d},\ldots,v_{d-1,d-1},v_{d-1,d},v_{dd}) ^{\top} \\
& := (\Sigma_{11},2\Sigma_{12},\ldots,2\Sigma_{1d},\Sigma_{22},2\Sigma_{23},\ldots,2\Sigma_{2d},\ldots,\Sigma_{d-1,d-1},2\Sigma_{d-1,d},\Sigma_{dd}) ^{\top}.
\end{aligned}
\end{equation*}
Define the density of $N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ parameterized in terms of $\boldsymbol{\mu}$ and $\boldsymbol{v}$ as
\begin{equation} \label{fv_defn}
f_v(\boldsymbol{\mu},\boldsymbol{v}) : = f(\boldsymbol{\mu},\boldsymbol{S}(\boldsymbol{v})), \quad \text{ where }\ S_{ij}(\boldsymbol{v}):= \begin{cases}
v_{ii} & \text{if } i=j , \\
v_{ij}/2 & \text{if } i \neq j .
\end{cases}
\end{equation}
For a $d\times d$ symmetric matrix $\boldsymbol{A}$, define a function $\boldsymbol{w}(\boldsymbol{A}) \in \mathbb{R}^{d(d+1)/2}$ that collects the unique elements of $\boldsymbol{A}$ as
\begin{align*}
\boldsymbol{w}(\boldsymbol{A}) & := (A_{11},2A_{12},\ldots,2A_{1d},A_{22},2A_{23},\ldots,2A_{2d},\ldots,A_{d-1,d-1},2A_{d-1,d},A_{dd}) ^{\top}.
\end{align*}
Then $f_v(\boldsymbol{\mu},\boldsymbol{v})$ and $f(\boldsymbol{\mu},\boldsymbol{\Sigma})$ are related as
\[
f(\boldsymbol{\mu},\boldsymbol{\Sigma}) = f_v(\boldsymbol{\mu},\boldsymbol{w}(\boldsymbol{\Sigma})).
\]

We introduce the following one-to-one mapping between $(\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\boldsymbol{v}_1,\boldsymbol{v}_2)$ and the reparameterized parameter $(\boldsymbol{\lambda}_{\boldsymbol{\mu}},\boldsymbol{\nu}_{\boldsymbol\mu},\boldsymbol{\lambda}_{\boldsymbol{v}},\boldsymbol{\nu}_{\boldsymbol{v}})$:
\begin{equation} \label{repara2}
\begin{pmatrix}
{\boldsymbol{\mu}}_1\\
{\boldsymbol{\mu}}_2\\
\boldsymbol{v}_1\\
\boldsymbol{v}_2
\end{pmatrix}
 =
\begin{pmatrix}
\boldsymbol{\nu}_{\boldsymbol\mu} + (1-\alpha) \boldsymbol{\lambda}_{\boldsymbol\mu} \\
\boldsymbol{\nu}_{\boldsymbol\mu} -\alpha \boldsymbol{\lambda}_{\boldsymbol\mu}\\
\boldsymbol{\nu}_{\boldsymbol{v}} + (1- \alpha)(2\boldsymbol{\lambda}_{\boldsymbol{v}}+ C_1 \boldsymbol{w}(\boldsymbol{\lambda}_{\boldsymbol\mu}\boldsymbol{\lambda}_{\boldsymbol\mu}^{\top}) )\\
\boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(2\boldsymbol{\lambda}_{\boldsymbol{v}}+ C_2 \boldsymbol{w}(\boldsymbol{\lambda}_{\boldsymbol\mu}\boldsymbol{\lambda}_{\boldsymbol\mu}^{\top})
\end{pmatrix},
\end{equation}
where $C_1 := -(1/3)(1 + \alpha)$ and $C_2 := (1/3)(2 - \alpha)$. Collect the reparameterized parameters, except for $\alpha$, into one vector $\boldsymbol{\psi}$ defined as
\begin{equation*}
\boldsymbol{\psi}:=(\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu},\boldsymbol{\nu}_{\boldsymbol{v}},\boldsymbol{\lambda}_{\boldsymbol\mu},\boldsymbol{\lambda}_{\boldsymbol{v}}) \in \Theta_{\boldsymbol{\psi}}.
\end{equation*}
In the reparameterized model, the null hypothesis of $H_{01}:(\boldsymbol{\mu}_{1},\boldsymbol{v}_1) = (\boldsymbol{\mu}_{2},\boldsymbol{v}_2)$ is written as $H_{01}:(\boldsymbol{\lambda}_{\boldsymbol\mu},\boldsymbol{\lambda}_{\boldsymbol v})= \boldsymbol{0}$, and the density is given by
\begin{equation} \label{loglike}
\begin{aligned}
g(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\psi},\alpha) & = \alpha f_v\left(\boldsymbol{x}\middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu}+(1-\alpha)\boldsymbol{\lambda}_{\boldsymbol\mu}, \boldsymbol{\nu}_{\boldsymbol{v}} + (1 - \alpha)(2\boldsymbol{\lambda}_{\boldsymbol{v}} + C_1 \boldsymbol{w}(\boldsymbol{\lambda}_{\boldsymbol\mu}\boldsymbol{\lambda}_{\boldsymbol\mu}^{\top}) ) \right) \\
& \quad + (1 - \alpha) f_v \left(\boldsymbol{x} \middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu} -\alpha\boldsymbol{\lambda}_{\boldsymbol\mu},\boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(2\boldsymbol{\lambda}_{\boldsymbol v}+ C_2 \boldsymbol{w}( \boldsymbol{\lambda}_{\boldsymbol\mu}\boldsymbol{\lambda}_{\boldsymbol\mu}^{\top}) ) \right).
\end{aligned}
\end{equation}


Partition $\boldsymbol{\psi}$ as $\boldsymbol{\psi} = (\boldsymbol{\eta}^{\top},\boldsymbol{\lambda}^{\top})^{\top}$, where $\boldsymbol{\eta}: = (\boldsymbol{\gamma}^{\top},\boldsymbol{\nu}_{\boldsymbol\mu}^{\top},\boldsymbol{\nu}_{\boldsymbol{v}}^{\top})^{\top} \in \Theta_{\boldsymbol{\eta}}$ and $\boldsymbol{\lambda}: = (\boldsymbol{\lambda}_{\boldsymbol\mu}^{\top},\boldsymbol{\lambda}_{\boldsymbol v}^{\top})^{\top} \in\Theta_{\boldsymbol{\lambda}}$. Denote the true values of $\boldsymbol{\eta}$, $\boldsymbol{\lambda}$, and $\boldsymbol{\psi}$ by $\boldsymbol{\eta}^*: = ((\boldsymbol{\gamma}^*)^{\top},({\boldsymbol{\mu}}^*)^{\top},(\boldsymbol{v}^{*})^{\top})^{\top}$, $\boldsymbol{\lambda}^*: = \boldsymbol{0}$, and $\boldsymbol{\psi}^* = ((\boldsymbol{\eta}^*)^{\top}, \boldsymbol{0}^{\top})^{\top}$, respectively.  Under this reparameterization, the first derivative of (\ref{loglike}) with respect to (w.r.t., hereafter) $\boldsymbol{\eta}$ under $\boldsymbol{\psi} = \boldsymbol{\psi}^*$ is identical to the first derivative of the density of the one-component model:
\begin{equation}\label{dvareta}
\nabla_{\boldsymbol{\eta}}g(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\psi}^*,\alpha) =\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top},\boldsymbol{v}^{\top})^{\top}} f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}^*,\boldsymbol{v}^{*}).
\end{equation}
On the other hand,   the first, second, and third derivatives of $g(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\psi},\alpha)$ w.r.t.\ $\boldsymbol{\lambda}_{\boldsymbol\mu}$ and the first derivative w.r.t.\ $\boldsymbol{\lambda}_{\boldsymbol v}$ become zero when evaluated at $\boldsymbol{\psi}=\boldsymbol{\psi}^*$. Consequently, the information on $\boldsymbol{\lambda}_{\boldsymbol\mu}$ and $\boldsymbol{\lambda}_{\boldsymbol v}$ is provided by the fourth derivative w.r.t.\ $\boldsymbol{\lambda}_{\boldsymbol\mu}$, the cross-derivative w.r.t.\ $\boldsymbol{\lambda}_{\boldsymbol{\mu}}$ and $\boldsymbol{\lambda}_{\boldsymbol v}$, and the second derivative w.r.t.\ $\boldsymbol{\lambda}_{\boldsymbol v}$.

We derive the asymptotic distribution of the LRTS.  Let $f^*_v$ and $\nabla f^*_v$ denote $f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}^*,\boldsymbol{v}^*)$ and $\nabla f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}^*,\boldsymbol{v}^*)$, and let $d_\eta:=(p+d+d(d+1)/2)$, $d_{\mu v}:=d(d+1)(d+2)/6$, and $d_{\mu^4}:=d(d+1)(d+2)(d+3)/24$. Define the score vector $\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z})$ as
\begin{equation} \label{score_defn}
\begin{aligned}
\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z}) & :=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\eta}}\\
\boldsymbol{s}_{\boldsymbol{\lambda}}
\end{pmatrix} :=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\eta}}\\
\boldsymbol{s}_{\boldsymbol{\mu v}}\\
\boldsymbol{s}_{\boldsymbol{\mu}^4}
\end{pmatrix}
\quad \text{with }\
\underset{(d_\mu \times 1)}{\boldsymbol{s}_{\boldsymbol{\eta}}} := \frac{\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top},\boldsymbol{v}^{\top})^{\top}}f^*_v }{ f^*_v}, \\
\underset{(d_{\mu v} \times 1)}{\boldsymbol{s}_{\boldsymbol{\mu v}}} &:= \left\{\frac{\nabla_{\mu_i \mu_j \mu_k} f^*_v}{ 3! f^*_v} \right\}_{1 \leq i \leq j \leq k \leq d}, \quad
\underset{(d_{\mu^4} \times 1)}{\boldsymbol{s}_{\boldsymbol{\mu}^4}} := \left\{ \frac{\nabla_{\mu_i \mu_j \mu_k \mu_\ell} f^*_v }{ 4! f^*_v} \right\}_{1 \leq i \leq j \leq k \leq \ell \leq d} ,
\end{aligned}
\end{equation}
where we suppress the dependence of $(\boldsymbol{s}_{\boldsymbol{\eta}}, \boldsymbol{s}_{\boldsymbol{\mu v}}, \boldsymbol{s}_{\boldsymbol{\mu}^4})$ on $(\boldsymbol{x},\boldsymbol{z})$. Collect the relevant reparameterized parameters as
\begin{equation} \label{tpsi_defn}
\boldsymbol{t}(\boldsymbol{\psi},\alpha) :=
\begin{pmatrix}
\boldsymbol{\eta}-\boldsymbol{\eta}^*\\
\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)
\end{pmatrix}
:=
\begin{pmatrix}
\boldsymbol{\eta}-\boldsymbol{\eta}^*\\
\alpha(1-\alpha)12 \boldsymbol{\lambda}_{\boldsymbol{\mu v}}\\
\alpha(1-\alpha)[12\boldsymbol{\lambda}_{\boldsymbol{v}^2}+ b(\alpha) \boldsymbol{\lambda}_{\boldsymbol{\mu}^4}]
\end{pmatrix},
\end{equation}
with $b(\alpha): = -(2/3) (\alpha^2 - \alpha + 1)<0$ and
\begin{equation} \label{lambda_muv_defn}
\begin{aligned}
\underset{(d_{\mu v} \times 1)}{\boldsymbol{\lambda}_{\boldsymbol{\mu v}}} &:= \{(\boldsymbol{\lambda}_{\boldsymbol{\mu v}})_{ijk}\}_{1 \leq i \leq j \leq k \leq d}, \text{ where } (\boldsymbol{\lambda}_{\boldsymbol{\mu v}})_{ijk} :=\sum_{(t_1,t_2,t_3) \in p_{12}(i,j,k)} \lambda_{\mu_{t_1}}\lambda_{v_{t_2t_3}}, \\
\underset{(d_{\mu^4} \times 1)}{\boldsymbol{\lambda}_{\boldsymbol{v}^2}} &:= \{(\boldsymbol{\lambda}_{\boldsymbol{v}^2})_{ijk\ell}\}_{1 \leq i \leq j \leq k \leq \ell \leq d}, \text{ where } (\boldsymbol{\lambda}_{\boldsymbol{v}^2})_{ijk\ell}:= \sum_{(t_1,t_2,t_3,t_4) \in p_{22}(i,j,k,\ell)} \lambda_{v_{t_1t_2}}\lambda_{v_{t_3t_4}} ,\\
\underset{(d_{\mu^4} \times 1)}{\boldsymbol{\lambda}_{\boldsymbol{\mu}^4}} &:= \{(\boldsymbol{\lambda}_{\boldsymbol{\mu}^4})_{ijk\ell}\}_{1 \leq i \leq j \leq k \leq \ell \leq d}, \text{ where } (\boldsymbol{\lambda}_{\boldsymbol{\mu}^4})_{ijk\ell}:= \sum_{(t_1,t_2,t_3,t_4) \in p(i,j,k,\ell)} \lambda_{\mu_{t_1}}\lambda_{\mu_{t_2}}\lambda_{\mu_{t_3}}\lambda_{\mu_{t_4}} ,
\end{aligned}
\end{equation}
where $\sum_{(t_1,t_2,t_3) \in p_{12}(i,j,k)}$ denotes the sum over all distinct permutations of $(i,j,k)$ to $(t_1,t_2,t_3)$ with $t_2 \leq t_3$, $\sum_{(t_1,t_2,t_3,t_4) \in p_{22}(i,j,k,\ell)}$ denotes the sum over all distinct permutations of $(i,j,k,\ell)$ to $(t_1,t_2,t_3,t_4)$ with $t_1 \leq t_2$ and $t_3 \leq t_4$, and $\sum_{(t_1,t_2,t_3,t_4) \in p(i,j,k,\ell)}$ denotes the sum over all distinct permutations of $(i,j,k,\ell)$ to $(t_1,t_2,t_3,t_4)$.  In (\ref{lambda_muv_defn}), $\boldsymbol{\lambda}_{\boldsymbol{\mu v}}$ is a function of $\boldsymbol{\lambda}_{\boldsymbol{\mu}} \otimes \boldsymbol{\lambda}_{\boldsymbol{v}}$ and corresponds to the score vector $\boldsymbol{s}_{\boldsymbol{\mu v}}$. $\boldsymbol{\lambda}_{\boldsymbol{v}^2}$ is a function of $\boldsymbol{\lambda}_{\boldsymbol{v}}^{\otimes 2}$, and $\boldsymbol{\lambda}_{\boldsymbol{\mu}^4}$ depends on $\boldsymbol{\lambda}_{\boldsymbol{\mu}}^{\otimes 4}$. Here, $\alpha(1-\alpha)12\boldsymbol{\lambda}_{\boldsymbol{\mu v}}^{\top} \boldsymbol{s}_{\boldsymbol{\mu v}}$ collects the unique elements that correspond to the cross-derivative with respect to $\boldsymbol{\lambda}_{\boldsymbol{\mu}}$ and $\boldsymbol{\lambda}_{\boldsymbol{v}}$ in the expansion of the log-likelihood function, and $\alpha(1-\alpha)[12\boldsymbol{\lambda}_{\boldsymbol{v}^2}+ b(\alpha) \boldsymbol{\lambda}_{\boldsymbol{\mu}^4}]^{\top} \boldsymbol{s}_{\boldsymbol{\mu}^4}$ collects the unique elements of the second-order terms with respect to $\boldsymbol{\lambda}_{\boldsymbol{v}}$ and the fourth-order terms with respect to $\boldsymbol{\lambda}_{\boldsymbol{\mu}}$.

Let $L_n(\boldsymbol{\psi},\alpha): = \sum_{i = 1}^n \log g(\boldsymbol{X}_i|\boldsymbol{Z}_i;\boldsymbol{\psi},\alpha)$ denote the reparameterized log-likelihood function. Let $\widehat{\boldsymbol{\psi}} : = \arg\max_{\boldsymbol{\psi} \in \Theta_{\boldsymbol{\psi}}} PL_n(\boldsymbol{\psi},\alpha)$ denote the PMLE of $\boldsymbol{\psi}$, where $\Theta_{\boldsymbol{\psi}}$ is defined so that the value of ${\boldsymbol{\vartheta}}_2$ implied by $\boldsymbol{\psi}$ is in $\Theta_{{\boldsymbol{\vartheta}}_2}$. Let $(\widehat{\boldsymbol{\gamma}}_0,\widehat{\boldsymbol{\mu}}_0,\widehat{\boldsymbol{\Sigma}}_0)$ denote the one-component MLE that maximizes the one-component log-likelihood function $L_{0,n}(\boldsymbol{\gamma},{\boldsymbol{\mu}},\boldsymbol{\Sigma}) := \sum_{i=1}^n \log f (\boldsymbol{X}_i|\boldsymbol{Z}_i;\boldsymbol{\gamma},{\boldsymbol{\mu}},\boldsymbol{\Sigma})$. Define the LRTS for testing $H_{01}$ as, with $\epsilon_1 \in (0,1/2)$,
\begin{equation} \label{LRT1}
LR_{n}(\epsilon_1) := \max_{\alpha \in [\epsilon_1,1-\epsilon_1]} 2\{L_n(\widehat{\boldsymbol{\psi}},\alpha) - L_{0,n}(\widehat{\boldsymbol{\gamma}}_0,\widehat{\boldsymbol{\mu}}_0,\widehat{\boldsymbol{\Sigma}}_0)\}.
\end{equation}
We could use the penalized LRTS defined by $PLR_{n}(\epsilon_1):= \max_{\alpha \in [\epsilon_1,1-\epsilon_1]} 2\{PL_n(\widehat{\boldsymbol{\psi}},\alpha) - L_{0,n}(\widehat{\boldsymbol{\gamma}}_0,\widehat{\boldsymbol{\mu}}_0,\widehat{\boldsymbol{\Sigma}}_0)\}$ instead of $LR_n(\epsilon_1)$. Because the effect of the penalty term is negligible under our assumptions, $PLR_n(\epsilon_1)$ has the same asymptotic distribution as $LR_n(\epsilon_1)$.

With $(\boldsymbol{s}_{\boldsymbol{\eta}}, \boldsymbol{s}_{\boldsymbol{\lambda}})$ defined in (\ref{score_defn}), define
\begin{equation} \label{I_lambda}
\begin{aligned}
&\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}} := E[\boldsymbol{s}_{\boldsymbol{\eta}}\boldsymbol{s}_{\boldsymbol{\eta}}^{\top}], \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}} := E[\boldsymbol{s}_{\boldsymbol{\lambda}}\boldsymbol{s}_{\boldsymbol{\lambda}}^{\top}], \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda\eta} } := E[\boldsymbol{s}_{\boldsymbol{\lambda}}\boldsymbol{s}_{\boldsymbol{\eta}}^{\top}],\\
&\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta \lambda}} := \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda \eta}}^{\top}, \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}:=\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}}-\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda\eta}}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}}^{-1}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta\lambda}}, \quad \boldsymbol{Z}_{\boldsymbol{\lambda}}:=(\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})^{-1} \boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}},
\end{aligned}
\end{equation}
where $\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \sim N(0,\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})$. The following sets characterize the limit of possible values of $\sqrt{n}\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)$ defined in (\ref{tpsi_defn}) as $n\rightarrow\infty$. Define
\begin{equation} \label{Lambda-e}
\begin{aligned}
\Lambda_{\boldsymbol{\lambda} }^{1} & := \left\{ \left((\boldsymbol{t}_{\boldsymbol{\mu v}} )^{\top}, ( \boldsymbol{t}_{\boldsymbol{\mu}^4} ) ^{\top} \right)^{\top} \in \mathbb{R}^{d_{\mu v}+d_{\mu^4}}: \boldsymbol{t}_{\boldsymbol{\mu v}} = \boldsymbol{\lambda}_{\boldsymbol{\mu v}},\ \boldsymbol{t}_{\boldsymbol{\mu}^4}= \boldsymbol{\lambda}_{\boldsymbol{v}^2} \text{ for some $\boldsymbol{\lambda} \in \mathbb{R}^{d+d(d+1)/2}$} \right\},\\
\Lambda_{\boldsymbol{\lambda} }^{2} & := \left\{ \left((\boldsymbol{t}_{\boldsymbol{\mu v}} )^{\top}, ( \boldsymbol{t}_{\boldsymbol{\mu}^4} ) ^{\top} \right)^{\top} \in \mathbb{R}^{d_{\mu v}+d_{\mu^4}}: \boldsymbol{t}_{\boldsymbol{\mu v}} = \boldsymbol{\lambda}_{\boldsymbol{\mu v}},\ \boldsymbol{t}_{\boldsymbol{\mu}^4}= -\boldsymbol{\lambda}_{\boldsymbol{\mu}^4} \text{ for some $\boldsymbol{\lambda} \in \mathbb{R}^{d+d(d+1)/2}$} \right\}.
\end{aligned}
\end{equation}
For $j=1,2$, define $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}$ by
\begin{equation} \label{t-lambda}
r(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}) = \inf_{\boldsymbol{t}_{\boldsymbol{\lambda}} \in \Lambda_{\boldsymbol{\lambda} }^{j}}r(\boldsymbol{t}_{\boldsymbol{\lambda}}), \quad r(\boldsymbol{t}_{\boldsymbol{\lambda}}) := (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}),
\end{equation}
where $\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}$, $\boldsymbol{Z}_{\boldsymbol{\lambda}}$, and $\Lambda_{\boldsymbol{\lambda} }^{j}$ for $j=1,2$ are defined in (\ref{I_lambda})--(\ref{Lambda-e}).

The following proposition establishes the asymptotic null distribution of the LRTS.
\begin{assumption} \label{A-taylor1}
$\boldsymbol{Z}$ has finite tenth moment.
\end{assumption}
\begin{proposition} \label{P-LR-N1}
Suppose that Assumptions \ref{assn_consis} and \ref{A-taylor1} hold, $a_n$ in (\ref{pen_pmle}) satisfies $a_n = O(1)$, and $\boldsymbol{\mathcal{I}}:=E[\boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z})\boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z})^{\top}]$ is finite and nonsingular. Then, under the null hypothesis of $M=1$,  $LR_{n}(\epsilon_1) \rightarrow_d \max\left\{ (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{1})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{1}, (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{2})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{2} \right\}$, where $LR_{n}(\epsilon_1)$ and $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}$ are defined in (\ref{LRT1}) and (\ref{t-lambda}), respectively.
\end{proposition}
For each $j=1,2$, the random variable $(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}$ is a projection of a Gaussian random variable on a cone $\Lambda_{\boldsymbol{\lambda} }^{j}$.

\begin{example}
When $d=1$ with $\boldsymbol{\lambda}=(\lambda_\mu,\lambda_v)^{\top}$, we have $\Lambda_{\boldsymbol{\lambda} }^{1} \cup \Lambda_{\boldsymbol{\lambda} }^{2} = \mathbb{R}^2$ and  $LR_{n}(\epsilon_1) \rightarrow_d \boldsymbol{Z}_{\boldsymbol{\lambda}}^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \boldsymbol{Z}_{\boldsymbol{\lambda}} \sim \chi^2(2)$.
When $d=2$, we have $\boldsymbol{\lambda}=(\boldsymbol{\lambda}_{\boldsymbol{\mu}}^{\top},\boldsymbol{\lambda}_{\boldsymbol{v}}^{\top})^{\top}=(\lambda_{\mu_1},\lambda_{\mu_2},\lambda_{v_{11}},\lambda_{v_{12}}, \lambda_{v_{22}})^{\top}$,
\[
\begin{aligned}
\boldsymbol{s}_{\boldsymbol{\mu v}} &= (\nabla_{\mu_1^3} f^*_v,\nabla_{\mu_1^2\mu_2} f^*_v,\nabla_{\mu_1\mu_2^2} f^*_v,\nabla_{\mu_2^3} f^*_v)^{\top}/ 3! f^*_v,  \\
\boldsymbol{s}_{\boldsymbol{\mu}^4} & = (\nabla_{\mu_1^4} f^*_v,\nabla_{\mu_1^3\mu_2} f^*_v,\nabla_{\mu_1^2\mu_2^2} f^*_v,\nabla_{\mu_1\mu_2^3} f^*_v,\nabla_{\mu_2^4} f^*_v)^{\top}/ 4! f^*_v,
 \end{aligned}
 \]
 and  $\Lambda_{\boldsymbol{\lambda} }^{j}$ is given by (\ref{Lambda-e}) with
\[
\begin{aligned}
\boldsymbol{t}_{\boldsymbol{\mu v}} & =(\lambda_{\mu_1}\lambda_{v_{11}},\lambda_{\mu_1}\lambda_{v_{12}}+\lambda_{\mu_2}\lambda_{v_{11}}, \lambda_{\mu_1}\lambda_{v_{22}}+\lambda_{\mu_2}\lambda_{v_{12}},\lambda_{\mu_2}\lambda_{v_{22}})^{\top},\\
 \boldsymbol{t}_{\boldsymbol{\mu }^4} &=
 \begin{cases}
 ( \lambda_{v_{11}}^2, 2\lambda_{v_{11}}\lambda_{v_{12}}, 2\lambda_{v_{11}}\lambda_{v_{22}}+\lambda_{v_{12}}^2, 2\lambda_{v_{12}}\lambda_{v_{22}} , \lambda_{v_{22}}^2)^{\top}& \text{if } j=1,\\
 -(\lambda_{\mu_1}^4,4\lambda_{\mu_1}^3\lambda_{\mu_2},6\lambda_{\mu_1}^2\lambda_{\mu_2}^2,4\lambda_{\mu_1}\lambda_{\mu_2}^3,\lambda_{\mu_2}^4)^{\top}& \text{if } j=2.
 \end{cases}
 \end{aligned}
 \]
\end{example}


\subsection{Likelihood ratio test of $H_0: M = M_0$ against $H_A: M = M_0 + 1$ for $M_0\geq 2$} \label{sec-general}

This section establishes the asymptotic distribution of the LRTS for testing the null hypothesis of $M_0$ components against the alternative of $M_0+1$ components for general $M_0 \geq 1$.


We consider a random sample of $n$ independent observations $\{\boldsymbol{X}_i,\boldsymbol{Z}_i\}_{i = 1}^n$ from the $M_0$-component $d$-variate finite normal mixture distribution, whose density with the true parameter value ${\boldsymbol{\vartheta}}_{M_0}^*=(\alpha_{1}^*,\ldots,\alpha_{M_0 - 1}^*,{\boldsymbol{\gamma}}^*,\boldsymbol{\mu}_1^*,\ldots,\boldsymbol{\mu}_{M}^*,\boldsymbol{\Sigma}_{1}^{*},\ldots,\boldsymbol{\Sigma}_{M}^{*})$ is
\begin{equation}
f_{M_0}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0}^*):=\sum_{j=1}^{M_0} \alpha_{j}^{*} f(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\gamma}}^*, \boldsymbol{\mu}_{j}^{*},\boldsymbol{\Sigma}_{j}^{*}), \label{true_model}
\end{equation}
where $\alpha_j^*>0$. We assume $(\boldsymbol{\mu}_{1}^*,\boldsymbol{\Sigma}^{*}_{1}) <\ldots< (\boldsymbol{\mu}_{M_0}^*,\boldsymbol{\Sigma}^{*}_{M_0})$ for identification. Let the density of an $(M_0+1)$-component mixture model be
\begin{equation}
f_{M_0+1}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0+1}):=\sum_{j=1}^{M_0+1}\alpha_j f(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j),\label{fitted_model}
\end{equation}
where ${\boldsymbol{\vartheta}}_{M_0+1} = (\alpha_1,\ldots,\alpha_{M_0},\boldsymbol{\gamma},\boldsymbol{\mu}_1,\ldots.,\boldsymbol{\mu}_{M_0+1},\boldsymbol{\Sigma}_1,\ldots,\boldsymbol{\Sigma}_{M_0+1})$. As in the case of the test of homogeneity, we partition the null hypothesis into two as $H_0 = H_{01} \cup H_{02}$, where $H_{01}: = \cup_{m=1}^{M_0} H_{0,1m}$ and $H_{02} := \cup_{m=1}^{M_0+1} H_{0,2m}$ with
\[
H_{0,1m} :(\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1) < \cdots < (\boldsymbol{\mu}_{m},\boldsymbol{\Sigma}_{m}) = (\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1}) < \cdots < (\boldsymbol{\mu}_{M_0 + 1},\boldsymbol{\Sigma}_{M_0 + 1}) \ \text{and}\ H_{0,2m}: \alpha_m = 0.
\]
The inequality constraints are imposed on $(\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)$ for identification.

We focus on testing $H_{01}$ because the LRTS for testing $H_{02}$ has infinite Fisher information unless a stringent restriction is imposed on the admissible values of $\boldsymbol{\Sigma}_j$ \citep{kasaharashimotsu15jasa}. Define the set of values of ${\boldsymbol{\vartheta}}_{M_0+1}$ that yields the true density (\ref{true_model}) as
\begin{equation*}
\Upsilon^*:=\{{\boldsymbol{\vartheta}}_{M_0 + 1}: f_{M_0 + 1}(\boldsymbol{X}|\boldsymbol{Z};{\boldsymbol{\vartheta}}_{M_0 + 1}) = f_{M_0}(\boldsymbol{X}|\boldsymbol{Z};{\boldsymbol{\vartheta}}_{M_0}^*)\text{ with probability one}\}.
\end{equation*}
Under $H_{0,1m}$, the $(M_0 + 1)$-component model (\ref{fitted_model}) generates the true $M_0$-component density (\ref{true_model}) when $(\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m) = (\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1})=(\boldsymbol{\mu}_{m}^{*},\boldsymbol{\Sigma}^{*}_{m})$. Define the subset of $\Upsilon^*$ corresponding to $H_{0,1m}$ as
\begin{align}
\Upsilon_{1m}^* &:= \left\{{\boldsymbol{\vartheta}}_{M_0 + 1} \in \Theta_{{\boldsymbol{\vartheta}}_{M_0 + 1}}:\  \alpha_j = \alpha_{j}^{*}\ \text{and}\
(\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j)=(\boldsymbol{\mu}_{j}^*,\boldsymbol{\Sigma}_{j}^{*})\ \text{for $j < m$}; \right. \nonumber \\
& \qquad \left.
\alpha_m + \alpha_{m + 1} = \alpha_{m}^{*}\ \text{and}\ (\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)=(\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1})=(\boldsymbol{\mu}_{m}^{*},\boldsymbol{\Sigma}_{m}^{*});\   \right.\nonumber \\
& \qquad \left. \alpha_{j}=\alpha_{j - 1}^*\ \text{and}\ (\boldsymbol{\mu}_{j},\boldsymbol{\Sigma}_j)=(\boldsymbol{\mu}_{j - 1}^*,\boldsymbol{\Sigma}_{j - 1}^{*})\ \text{for $j> m + 1$};\ {\boldsymbol{\gamma}}={\boldsymbol{\gamma}}^* \right\}, \nonumber
\end{align}
and define $\Upsilon_1^*:= \Upsilon_{11}^* \cup \cdots \cup \Upsilon_{1M_0}^*$.

Let $\Theta_{{\boldsymbol{\vartheta}}_{M_0 + 1}}(\epsilon_1)$ be a subset of $\Theta_{{\boldsymbol{\vartheta}}_{M_0 + 1}}$ such that $\alpha_j\in[\epsilon_1,1-\epsilon_1]$ for $j = 1,\ldots,M_0 + 1$, and define the PMLE by
\begin{equation}\label{PMLEs}
\begin{aligned}
\widehat{\boldsymbol{\vartheta}}_{M_0+1}(\epsilon_1) & := \mathop{\arg \max}_{{\boldsymbol{\vartheta}}_{M_0+1}\in\Theta_{{\boldsymbol{\vartheta}}_{M_0+1}}(\epsilon_1)}PL_n({\boldsymbol{\vartheta}}_{M_0 + 1}),\\
\widehat{\boldsymbol{\vartheta}}_{M_0} & := \mathop{\arg \max}_{{\boldsymbol{\vartheta}}_{M_0}\in\Theta_{{\boldsymbol{\vartheta}}_{M_0}}}PL_{0,n}({\boldsymbol{\vartheta}}_{M_0}),
\end{aligned}
\end{equation}
where $PL_n({\boldsymbol{\vartheta}}_{M_0 + 1}):=L_n({\boldsymbol{\vartheta}}_{M_0 + 1})+p_n(\boldsymbol{\vartheta}_{M_0+1})$ and $PL_{0,n}({\boldsymbol{\vartheta}}_{M_0}):=L_{0,n}({\boldsymbol{\vartheta}}_{M_0})+p_n(\boldsymbol{\vartheta}_{M_0})$ with $L_n({\boldsymbol{\vartheta}}_{M_0 + 1}):=\sum_{i = 1}^n \log f_{M_0 + 1}(\boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_{M_0 + 1})$ and $L_{0,n}({\boldsymbol{\vartheta}}_{M_0}):=\sum_{i = 1}^n \log f_{M_0}(\boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_{M_0})$ for the density (\ref{true_model})--(\ref{fitted_model}) and the penalty function in (\ref{pen_pmle}).
We consider the LRTS for testing $H_{01}$ given by
\begin{equation}\label{LR-M_0}
LR_{n}^{M_0}(\epsilon_1):= 2\{L_n(\widehat{{\boldsymbol{\vartheta}}}_{M_0+1}(\epsilon_1))- L_{0,n}(\widehat{{\boldsymbol{\vartheta}}}_{M_0})\}.
\end{equation}

Collect the score vector for testing $H_{0,11},\ldots,H_{0,1M_0}$ into one vector as
\begin{equation}
\widetilde{\boldsymbol{s}}(\boldsymbol{x},\boldsymbol{z}) :=
\begin{pmatrix}
\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}} \\ \widetilde{\boldsymbol{s}}_{\boldsymbol{\lambda}}
\end{pmatrix},\ \text{ where }\
\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}}:=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\alpha}} \\ \boldsymbol{s}_{(\boldsymbol{\gamma},\boldsymbol{\mu}, \boldsymbol{v})}
\end{pmatrix}\ \text{and}\ \widetilde{\boldsymbol{s}}_{\boldsymbol{\lambda}}: =
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\mu}\boldsymbol{v}}^1 \\
\boldsymbol{s}_{\boldsymbol{\mu}^4}^1 \\
 \vdots \\
\boldsymbol{s}_{\boldsymbol{\mu}\boldsymbol{v}}^{M_0} \\
\boldsymbol{s}_{\boldsymbol{\mu}^4}^{M_0}
\end{pmatrix}, \label{stilde}
\end{equation}
where, with $f_0^*:=f_{M_0}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0}^*)$ and for $m=1,\ldots,M_0$,
\begin{equation} \label{sh}
\begin{aligned}
\boldsymbol{s}_{\boldsymbol{\alpha}} & :=
\begin{pmatrix}
f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_1^{*},\boldsymbol{v}_1^*)-f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_{M_0}^{*},\boldsymbol{v}_{M_0}^*) \\
\vdots \\
f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_{M_0-1}^{*},\boldsymbol{v}_{M_0-1}^*)-f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}^{*}_{M_0},\boldsymbol{v}_{M_0}^*) \\
\end{pmatrix} \Bigl/ f_0^*, \bigr.
\\
\boldsymbol{s}_{(\boldsymbol{\gamma},\boldsymbol{\mu}, \boldsymbol{v})} & :=\sum_{m=1}^{M_0}\alpha_{m}^* \nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top}, \boldsymbol{v}^{\top})^{\top}} f_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_m^{*},\boldsymbol{v}_m^*)/ f_0^*, \\
\boldsymbol{s}_{\boldsymbol{\mu v}}^m &:= \left\{\alpha_{m}^* \nabla_{\mu_i \mu_j \mu_k} f^*_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_m^{*},\boldsymbol{v}_m^*) / (3! f_0^*) \right\}_{1 \leq i \leq j \leq k \leq d}, \\
\boldsymbol{s}_{\boldsymbol{\mu}^4}^m &:= \left\{ \alpha_{m}^* \nabla_{\mu_i \mu_j \mu_k \mu_\ell} f^*_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_m^{*},\boldsymbol{v}_m^*) / (4! f_0^*) \right\}_{1 \leq i \leq j \leq k \leq \ell \leq d} .
\end{aligned}
\end{equation}

Define
\begin{equation}\label{Itilde}
\begin{aligned}
&\widetilde{\boldsymbol{\mathcal{I}}}:= E[ \widetilde{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z})\widetilde{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z})^{\top}],\ \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}}:= E[\widetilde {\boldsymbol{s}}_{\boldsymbol{\eta}} \widetilde {\boldsymbol{s}}_{\boldsymbol{\eta}}^{\top}], \ \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}:= E[\widetilde {\boldsymbol{s}}_{\boldsymbol{\lambda}} \widetilde {\boldsymbol{s}}_{\boldsymbol{\eta}}^{\top}], \ \\
&\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta\lambda}}:= \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}^{\top},\quad \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}}:= E[\widetilde {\boldsymbol{s}}_{\boldsymbol{\lambda}} \widetilde {\boldsymbol{s}}_{\boldsymbol{\lambda}}^{\top}],\ \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} :=\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}} - \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}}^{-1} \widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta\lambda}}.
\end{aligned}
\end{equation}
Let $\widetilde{\boldsymbol{G}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}=((\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^{1})^{\top},\ldots,(\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^{M_0})^\top)^\top \sim N(0,\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})$ be an $\mathbb{R}^{M_0 (d_{\mu v} + d_{\mu^4})}$--valued random vector, and define ${\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m:=E[\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m (\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m)^\top]$ and $\boldsymbol{Z}_{\boldsymbol{\lambda}}^m:=({\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m)^{-1}\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m$. For $j=1,2$, similar to $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}$ in the test of homogeneity, define $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^{j}$ by
\begin{equation*}
r^m(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^{j}) = \inf_{\boldsymbol{t}_{\boldsymbol{\lambda}} \in \Lambda_{\boldsymbol{\lambda} }^{j}}r^m(\boldsymbol{t}_{\boldsymbol{\lambda}}), \quad r^m(\boldsymbol{t}_{\boldsymbol{\lambda}}) := (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}^m)^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}^m),
\end{equation*}
where $\Lambda_{\boldsymbol{\lambda} }^{j}$ is defined in (\ref{Lambda-e}).
The following proposition gives the asymptotic null distribution of the LRTS for testing $H_{01}$. In the neighborhood of $\Upsilon_{1h}^*$, the log-likelihood function permits a quadratic approximation in terms of polynomials of the parameters similar to testing $H_0:M=1$ against $H_A:M=2$. Consequently, the LRTS is asymptotically distributed as the maximum of $M_0$ random variables.

\begin{assumption}\label{A-vec-2} (a) $\alpha_{j}^*\in [\epsilon_1,1-\epsilon_1]$ for $j = 1,\ldots,M_0$. (b) $\widetilde{\boldsymbol{\mathcal{I}}}$ defined in (\ref{Itilde}) is nonsingular.
\end{assumption}

\begin{proposition} \label{local_lr-2}
Suppose that Assumptions \ref{assn_consis}, \ref{A-taylor1}, and \ref{A-vec-2} hold  and $a_n$ in (\ref{pen_pmle}) satisfies $a_n = o(1)$. Then, under the null hypothesis $H_0: M=M_0$, $LR_{n}^{M_0}(\epsilon_1) \rightarrow_d \max\{v_1,\ldots, v_{M_0}\}$, where $v_m := \max\left\{ (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^{1})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^1, (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^{2})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^2\right\}$.
\end{proposition}

\section{EM test} \label{section:emtest}

Implementing the likelihood ratio test in Section \ref{sec:hetero} requires the researcher to choose a lower bound $\epsilon_1$ on $\alpha_j$ and assume $\alpha_j^* >\epsilon_1$.  In this section, we develop an EM test of $H_0:M= M_0$ against $H_A:M = M_0 + 1$  that does not require such a lower bound on $\alpha_j$.  For brevity, we  suppress covariate $\boldsymbol{Z}$ in this section. First, we develop an EM test statistic for testing $H_{0,1m}: (\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)=(\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1})$. We construct $M_0$ sets $\{D_1^*,\cdots,D_{M_0}^*\}$ of admissible values of $(\boldsymbol{\mu},\boldsymbol{\Sigma})$, such that $D_m$ contains $(\boldsymbol{\mu}_m^*,\boldsymbol{\Sigma}_m^{*})$ but no other $(\boldsymbol{\mu}_j^*,\boldsymbol{\Sigma}_j^{*})$'s for $j \neq m$. For example, as in our simulation, we may assume that the first element of $\boldsymbol{\mu}$ are distinct, let $\overline{\mu}_j^*:=(\mu_{j1}^* +\mu_{j + 1,1}^*)/2$ with $\mu_{j1}$ denoting the first element of $\boldsymbol{\mu}_j$, and set $D_1^* = (-\infty , \overline{\mu}_{1}^*] \times\Theta_{\tilde{\boldsymbol{\mu}}}\times\Theta_{\boldsymbol{\Sigma}}$, $D_j^* =[\overline{\mu}_{j-1}^*, \overline{\mu}_{j}^*]\times\Theta_{\tilde{\boldsymbol{\mu}}}\times\Theta_{\boldsymbol{\Sigma}}$ for $j = 2,\ldots,M_0 - 1$, and $D_{M_0}^* = [\overline{\mu}_{M_0 - 1}^*, \infty ) \times\Theta_{\tilde{\boldsymbol{\mu}}} \times \Theta_{\boldsymbol{\Sigma}}$, where $\Theta_{\tilde{\boldsymbol{\mu}}}$ denotes the space of $\tilde{\boldsymbol{\mu}}:=(\mu_2,\ldots,\mu_d)\top$.

Collect the mixing parameters of the $(M_0 + 1)$-component model into one vector as $\boldsymbol{\varsigma}:=(\boldsymbol{\mu}_1,\ldots,\boldsymbol{\mu}_{M_0 + 1},\boldsymbol{\Sigma}_1,\ldots,\boldsymbol{\Sigma}_{M_0 + 1}) \in \Theta_{\boldsymbol{\varsigma}}:=\Theta_{\boldsymbol{\mu}}^{M_0 + 1}\times \Theta_{\boldsymbol{\Sigma}}^{M_0 + 1}$. For $m = 1,\ldots,M_0$, define a restricted parameter space of $\boldsymbol{\varsigma}$ by $\Xi_m^*: = \{\boldsymbol{\varsigma} \in \Theta_{\boldsymbol{\varsigma}}: (\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j) \in D_j^* \ \text{for}\ j = 1,\ldots, m - 1;\ (\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m),\ (\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1})\in D_m^*;\ (\boldsymbol{\mu}_j,\boldsymbol{\Sigma}_j) \in D_{j - 1}^* \ \text{for}\ j = m + 2,\ldots,M_0 + 1 \}$. Let $\widehat{\Xi}_m$ and $\widehat{D}_m$ be consistent estimates of $\Xi_m^*$ and $D_m^*$, which can be constructed from  the PMLE  of the $M_0$-component model. We test $H_{0,1m}: (\boldsymbol{\mu}_m,\boldsymbol{\Sigma}_m)=(\boldsymbol{\mu}_{m + 1},\boldsymbol{\Sigma}_{m + 1})$ by estimating the $(M_0 + 1)$-component model (\ref{fitted_model}) under the restriction $\boldsymbol{\varsigma} \in \widehat{\Xi}_m$. For example, when we test $H_{0,11}:(\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1) = (\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2)$ in a three-component model, the restriction can be given as $(\boldsymbol{\mu}_1,\boldsymbol{\Sigma}_1), (\boldsymbol{\mu}_2,\boldsymbol{\Sigma}_2) \in \widehat{D}_1$ and $(\boldsymbol{\mu}_3,\boldsymbol{\Sigma}_3) \in \widehat{D}_2$.

Define the penalty term $p_n^m(\boldsymbol{\vartheta}_{M_0+1})$ on $\boldsymbol{\Sigma}_j$'s as
\begin{equation}\label{penalty-em}
p_n^m(\boldsymbol{\vartheta}_{M_0+1}) := \sum_{j=1}^{M_0+1}p^m_{n}(\boldsymbol{\Sigma}_j;\boldsymbol{\Omega}_j), \quad p^m_{n}(\boldsymbol{\Sigma}_j;\boldsymbol{\Omega}_j):=- a_n \left\{ \text{tr}(\boldsymbol{\Omega}_j\boldsymbol{\Sigma}_j^{-1}) - \log ( \det({\boldsymbol{\Omega}}_j\boldsymbol{\Sigma}_j^{-1})) -d \right\},
\end{equation}
with $\boldsymbol{\Omega}_j = \widehat{\boldsymbol{\Sigma}}_j$ for $j=1,\ldots,m-1$, $\boldsymbol{\Omega}_j =\widehat{\boldsymbol{\Sigma}}_m$ for $j=m,m+1$, and $\boldsymbol{\Omega}_j =\widehat{\boldsymbol{\Sigma}}_{j-1}$ for $j=m+2,\ldots,M_0+1$,
where $\widehat{\boldsymbol{\Sigma}}_j$ is a consistent estimator of $\boldsymbol{\Sigma}_j$ from the $M_0$-component PMLE. This penalty term is a multivariate version of the one in \citet{chenlifu12jasa} and satisfies $p^m_{n}(\boldsymbol{\Omega}_j;\boldsymbol{\Omega}_j)=0$. Let $\mathcal{T}$ be a finite set of numbers from $(0,0.5]$, and let $p(\tau) \leq 0$ be a penalty term that is continuous in $\tau$, $p(0.5)=0$, and $p(\tau) \to -\infty$ as $\tau$ goes to $0$.

For each $\tau_0 \in \mathcal{T}$, define the restricted penalized MLE as $\boldsymbol{\vartheta}_{M_0+1}^{m(1)}(\tau_0) := \\\mathop{\arg \max}_{\boldsymbol{\vartheta}_{M_0+1} \in \Theta^m(\tau_0)} (PL_n^m({\boldsymbol{\vartheta}}_{M_0 + 1}) + p(\tau_0))$, where $\Theta^m(\tau) := \{ \boldsymbol{\vartheta}_{M_0+1} \in \Theta_{\boldsymbol{\vartheta}_{M_0+1}}: \alpha_{m}/(\alpha_{m} + \alpha_{m + 1})=\tau \text{ and } \boldsymbol{\varsigma}\in \hat \Xi_m\}$ and $PL_{n}^m({\boldsymbol{\vartheta}}_{M_0+1}):= \sum_{i=1}^n f_{M_0+1}(\boldsymbol{X}_i ;\boldsymbol{\vartheta}_{M_0+1}) + p_n^m(\boldsymbol{\vartheta}_{M_0+1})$. Starting from $\boldsymbol{\vartheta}_{M_0+1}^{m(1)}(\tau_0)$, we update $\boldsymbol{\vartheta}_{M_0+1}$ by the following generalized EM algorithm. Henceforth, we suppress $(\tau_0)$ from $\boldsymbol{\vartheta}_{M_0+1}^{m(k)}(\tau_0)$. Suppose we have already calculated $\boldsymbol{\vartheta}_{M_0+1}^{m(k)}$. For $i = 1,\ldots,n$ and $j = 1,\ldots,M_0+1$, define the weights for an E-step as
\begin{align*}
w_{ij}^{(k)} &:=
\begin{cases}
\alpha_j^{(k)} f(\boldsymbol{X}_i;{\boldsymbol{\mu}}_j^{(k)},\boldsymbol{\Sigma}_j^{(k)}) / f_{M_0+1}(\boldsymbol{X}_i;{\boldsymbol{\vartheta}}_{M_0+1}^{m(k)}) & \mbox{for } j=1,\ldots,m-1,\\
\alpha_{j}^{(k)} f(\boldsymbol{X}_i;{\boldsymbol{\mu}}_j^{(k)},\boldsymbol{\Sigma}_j^{(k)}) / f_{M_0+1}(\boldsymbol{X}_i;{\boldsymbol{\vartheta}}_{M_0+1}^{m(k)}) & \mbox{for } j=m+2,\ldots,M_0+1,\\
\end{cases} \\
w_{im}^{(k)} &:= \frac{\tau^{(k)}(\alpha_m^{(k)}+\alpha_{m+1}^{(k)}) f(\boldsymbol{X}_i;{\boldsymbol{\mu}}_m^{(k)},\boldsymbol{\Sigma}_m^{(k)})}{f_{M_0+1}(\boldsymbol{X}_i;{\boldsymbol{\vartheta}}_{M_0+1}^{m(k)})}, \ w_{i,m+1}^{(k)} := \frac{(1-\tau^{(k)})(\alpha_m^{(k)}+\alpha_{m+1}^{(k)})f(\boldsymbol{X}_i;{\boldsymbol{\mu}}_{m+1}^{(k)},\boldsymbol{\Sigma}_{m+1}^{(k)}) }{f_{M_0+1}(\boldsymbol{X}_i;{\boldsymbol{\vartheta}}_{M_0+1}^{m(k)})}.
\end{align*}
In an M-step, update $\tau$ and $\boldsymbol{\alpha}$ by
\begin{align*}
\tau^{(k+1)} &:= \mathop{\arg \max}_{\tau} \left\{\sum_{i=1}^n w_{im}^{(k)} \log(\tau) + \sum_{i=1}^n w_{i,m+1}^{(k)} \log(1-\tau) + p(\tau) \right\},\\
\alpha_j^{(k+1)} &:=
n^{-1} \sum_{i=1}^n w_{ij}^{(k)} \quad \mbox{for } j=1, \ldots,M_0+1,
\end{align*}
and update $\boldsymbol{\mu}_j$ and $\boldsymbol{\Sigma}_j$ for $j=1,\ldots,M_0+1$ by
\begin{equation*}
\begin{aligned}
\boldsymbol{\mu}_j^{(k+1)} & : = \frac{\sum_{i=1}^n w_{ij}^{(k)} \boldsymbol{X}_i}{\sum_{i=1}^n w_{ij}^{(k)} },\quad
\boldsymbol{\Sigma}_j^{(k+1)} : = \frac{2a_n \boldsymbol{\Omega}_j + \boldsymbol{S}_j^{(k+1)}}{2a_n + \sum_{i=1}^n w_{ij}^{(k)}},
 \\
\end{aligned}
\end{equation*}
where $\boldsymbol{S}_j^{(k+1)}:= \sum_{i=1}^nw_{ij}^{(k)} \left(\boldsymbol{X}_i-\boldsymbol{\mu}_j^{(k+1)}\right) \left(\boldsymbol{X}_i-\boldsymbol{\mu}_j^{(k+1)}\right)^{\top}$.
The penalized likelihood value never decreases after each generalized EM step \citep[Theorem 1]{dempster77jrssb}. Note that $\boldsymbol{\vartheta}_{M_0+1}^{m(k)}$ for $k \geq 2$ does not use the restriction $\hat \Xi_m$. For each $\tau_0 \in \mathcal{T}$ and $k$, define
\begin{equation*}
\text{M}_{n}^{m(k)}(\tau_0) : = 2\left\{PL_{n}^m(\boldsymbol{\vartheta}_{M_0+1}^{m(k)}(\tau_0) ) +p(\tau^{(k)}) - L_{0,n}(\widehat{\boldsymbol{\vartheta}}_{M_0}) \right\},
\end{equation*}
where $\widehat{\boldsymbol{\vartheta}}_{M_0}$ and $L_{0,n}({\boldsymbol{\vartheta}}_{M_0})$ are defined in (\ref{PMLEs}).

Finally, with a pre-specified number $K$, define the \textit{local EM test statistic} for testing $H_{0,1m}$ by taking the maximum of $\text{M}_n^{m(K)}(\tau_0)$ over $\tau_0\in\mathcal{T}$ as $\text{EM}_{n}^{m(K)} : = \max \{\text{M}_{n}^{m(K)}(\tau_0): \tau_0 \in \mathcal{T} \}$. The \textit{EM test statistic} is defined as the maximum of $M_0$ local EM test statistics:
\begin{equation}\label{EM-test}
\text{EM}_{n}^{(K)}: = \max\left\{\text{EM}_{n}^{1(K)},\text{EM}_{n}^{2(K)},\ldots,\text{EM}_{n}^{M_0(K)}\right\}.
\end{equation}
The following proposition shows that for any finite $K$, the EM test statistic is asymptotically equivalent to the penalized LRTS for testing $H_{01}$.


\begin{proposition} \label{EM_stat-1}
Suppose that  Assumptions \ref{assn_consis} and \ref{A-vec-2} hold, $a_n$ in (\ref{penalty-em}) satisfies $a_n = O(1)$, and $\{0.5\} \in \mathcal{T}$. Then, under the null hypothesis $H_0: M=M_0$, for any fixed finite $K$, $\text{EM}_{n}^{{(K)}} \rightarrow_d \max\{v_1,\ldots, v_{M_0}\}$ as $n \rightarrow \infty$, where the $v_m$'s are given in Proposition \ref{local_lr-2}.
\end{proposition}


\section{Asymptotic distribution under local alternatives}
In this section, we derive the asymptotic distribution of our LRTS  and EM test statistic  under local alternatives. For brevity, we focus on the case of testing $H_0: M=1$ against $H_A: M=2$.

Given a local parameter $\boldsymbol{h} =(\boldsymbol{h}_{\boldsymbol \eta}^{\top},
\boldsymbol{h}_{\boldsymbol \lambda }^{\top})^{\top}$ and $\alpha  \in (\epsilon_1,1-\epsilon_1)$, we consider the sequence of contiguous local alternatives $\boldsymbol{\vartheta}_{n} = (\boldsymbol{\psi}_n^{\top},\alpha_n)^{\top} = (\boldsymbol{\eta}_n^{\top},\boldsymbol\lambda_n^{\top},\alpha_n)^{\top}\in\Theta_{\boldsymbol{\eta}}\times\Theta_{\boldsymbol\lambda}\times\Theta_{\alpha}$ such that, with $\boldsymbol{t}_{\boldsymbol\lambda}(\boldsymbol\lambda,\alpha)$ given by (\ref{tpsi_defn}),
\begin{equation}\label{local-alternative}
\boldsymbol{h}_{\boldsymbol \eta}= \sqrt{n}(\boldsymbol{\eta}_n-\boldsymbol\eta^*),\quad \boldsymbol{h}_{\boldsymbol\lambda} =\sqrt{n} \boldsymbol{t}_{\boldsymbol\lambda}(\boldsymbol\lambda_{n},\alpha_n)+o(1),\quad \text{and }\ \alpha_n = \alpha+o(1).
\end{equation}

Let $\mathbb{P}_{\boldsymbol{\vartheta}}^n$ be the probability measure on $\{\boldsymbol{X}_i\}_{i=1}^n$ conditional on $\{\boldsymbol{Z}_i\}_{i=1}^n$ under $\boldsymbol{\vartheta}$. Then, for the density (\ref{loglike}), the log-likelihood ratio is given by
\[
\log \frac{d\mathbb{P}_{\boldsymbol{\vartheta}_n}^n}{d \mathbb{P}_{\boldsymbol{\vartheta}^*}^n} = L_n(\boldsymbol{\psi}_n,\alpha_n)-L_n(\psi^*,\alpha)= \sum_{i=1}^n \log \left( \frac{g(X_i|Z_i;\boldsymbol{\eta}_n,\boldsymbol{\lambda}_n,\alpha_n)}{g(X_i|Z_i;\boldsymbol{\eta}^*,\boldsymbol{0},\alpha)} \right).
\]

The following proposition provides the asymptotic distribution of the LRTS under contiguous local alternatives.
\begin{proposition} \label{P-LAN2} Suppose that the assumptions of Proposition \ref{P-LR-N1} hold. Consider a sequence of contiguous local alternatives $\boldsymbol{\vartheta}_{n} = ((\boldsymbol{\eta}^*)^{\top},(\boldsymbol{\lambda}_n) ^{\top},\alpha_n)$, where $\boldsymbol{\lambda}_n$ and $\alpha_n$  are given by (\ref{local-alternative}). Then, under $H_{1n}: \boldsymbol{\vartheta} = \boldsymbol{\vartheta}_{n}$, we have
$LR_{n}(\epsilon_1) \rightarrow_d \max\left\{ (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{1})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{1}, (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{2})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{2} \right\}$, where $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^1$ and $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^2$ are defined as in (\ref{t-lambda}) but replacing $\boldsymbol{Z_{\lambda}}$ with $\left(\boldsymbol{I_{\lambda,\eta}}\right)^{-1} \boldsymbol{G_{\lambda.\eta}}+\boldsymbol{h}_{\boldsymbol\lambda}$.
\end{proposition}

In this proposition, the local alternatives are implicitly defined through the condition that $\boldsymbol{h}_{\boldsymbol\lambda} =\sqrt{n} \boldsymbol{t}_{\boldsymbol\lambda} (\boldsymbol\lambda_{n} ,\alpha_n)+o(1)$. We now give an example for $d=1$, where we explicitly construct local alternatives with different orders, including the ones with order $n^{1/8}$.
\begin{example} When $d=1$, for $\overline\alpha\in(\epsilon_1,1-\epsilon_1)$ and $\overline{\boldsymbol{\lambda}}:=(\overline{{\lambda}}_{\mu},\overline{{\lambda}}_{v})^{\top}$, define
\begin{align*}
&H_{1n}^a : \boldsymbol{\vartheta}_n^a=((\boldsymbol{\eta}_n^a)^{\top},(\boldsymbol{\lambda}_{n}^a)^{\top},\alpha_n^a)^{\top} : = ((\boldsymbol\eta^*)^{\top},(\overline{\boldsymbol{\lambda}}/n^{1/4})^{\top},\overline\alpha + o(1))^{\top},\\
&H_{1n}^b : \boldsymbol{\vartheta}_n^b=((\boldsymbol{\eta}_n^b)^{\top},(\boldsymbol{\lambda}_{n}^b)^{\top},\alpha_n^b)^{\top} : = ((\boldsymbol\eta^*)^{\top},\overline{{\lambda}}_{\mu}/n^{1/8},\overline{{\lambda}}_{v}/n^{3/8},\overline\alpha + o(1))^{\top}.
\end{align*}
Then, for $j\in\{a,b\}$, $\boldsymbol{h}_{\boldsymbol\lambda}^j =\sqrt{n} \boldsymbol{t}_{\boldsymbol\lambda} (\boldsymbol\lambda_{n}^j ,\alpha_n^j)+o(1)$ holds with $\boldsymbol{h}_{\boldsymbol\lambda}^a := 12\overline\alpha(1-\overline\alpha)\times ( \overline{{\lambda}}_{\mu} \overline{{\lambda}}_{v}, \overline{{\lambda}}_{v}^2)$ and $\boldsymbol{h}_{\boldsymbol\lambda}^b := \overline\alpha(1-\overline\alpha)\times ( 12\overline{{\lambda}}_{\mu} \overline{{\lambda}}_{v}, b(\overline\alpha)\overline{{\lambda}}_{\mu}^4 )^{\top}$. Therefore, Proposition \ref{P-LAN2} gives the asymptotic distribution of $LR_{n}(\epsilon_1)$.
\end{example}


\section{Parametric bootstrap}

Given that it may not be easy to simulate the asymptotic distributions of the LRTS and the EM test statistic for testing $H_0: M=M_0$ against $H_A: M=M_0+1$, we provide the validity of parametric bootstrap. We consider the following parametric bootstrap to obtain the bootstrap critical value $c_{\alpha,B}$ and bootstrap $p$-value.
\begin{enumerate}
\item Using the observed data, compute $\widehat{\boldsymbol{\vartheta}}_{M_0}$ and compute $LR_{n}^{M_0}(\epsilon_1)$ in (\ref{LR-M_0}) and $\text{EM}_{n}^{{(K)}}$ in (\ref{EM-test}).
\item Given $\widehat{\boldsymbol{\vartheta}}_{M_0}$, generate $B$ independent samples $\{\boldsymbol X_1^b,\ldots,\boldsymbol X_n^b\}_{b=1}^B$ under $H_0$ with $\boldsymbol\vartheta_{M_0}=\widehat{\boldsymbol{\vartheta}}_{M_0}$ conditional on the observed value of $\{\boldsymbol Z_1,\ldots,\boldsymbol Z_n\}$.
\item For each simulated sample $\{\boldsymbol X_1^b,\ldots,\boldsymbol X_n^b\}$ with $\{\boldsymbol Z_1,\ldots,\boldsymbol Z_n\}$, compute $LR_{n}^{M_0,b}(\epsilon_1)$ and $\text{EM}_{n}^{{(K),b}}$ as in Step 1 for $b=1,\ldots,B$.
\item Let $c_{\alpha,B}$ be the $(1-\alpha)$ quantile of $\{LR_{n}^{M_0,b}\}_{b=1}^B$ or $\{EM_{n}^{(K),b}\}_{b=1}^B$, and define the bootstrap $p$-value as $B^{-1}\sum_{b=1}^B \mathbb{I}\{ LR_{n}^{M_0,b} > LR_{ n}^{M_0}\}$ or $B^{-1}\sum_{b=1}^B \mathbb{I}\{ EM_{n}^{(K),b} > EM_{ n}^{(K)}\}$.
\end{enumerate}
The following proposition shows the consistency of the bootstrap critical values $c_{\alpha,B}$ for testing $H_0: M=1$. The case of testing $H_0:M=M_0$ for $M_0\geq 2$ can be proven similarly.
\begin{proposition} \label{P-bootstrap}
Suppose that the assumptions of Proposition \ref{P-LR-N1} holds. Then, the bootstrap critical values $c_{\alpha,B}$ converge to the asymptotic critical values in probability as $n$ and $B$ go to infinity under $H_0$ and under the local alternatives described in Propositions \ref{P-LAN2}.
\end{proposition}


\section{  Homoscedastic multivariate finite normal mixture models}

In this section, we consider testing the order of homoscedastic multivariate normal mixtures. We consider the likelihood ratio test but do not consider the EM test because, unlike heteroscedastic normal mixtures, homoscedastic normal mixture models do not suffer from infinite Fisher information and unbounded likelihood.

\subsection{Likelihood ratio test of $H_0: M = 1$ against $H_A: M = 2$}

Consider a two-component normal mixture density function with common variance:
\begin{equation*}
f_2(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\vartheta}_2) := \alpha f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma} ,\boldsymbol{\mu}_1,\boldsymbol{\Sigma}) + (1-\alpha) f(\boldsymbol{x}|\boldsymbol{z}; \boldsymbol{\gamma},\boldsymbol{\mu}_2,\boldsymbol{\Sigma}),
\end{equation*}
with ${\boldsymbol{\vartheta}}_2 = (\alpha,\boldsymbol{\gamma}, \boldsymbol{\mu}_1, \boldsymbol{\mu}_2,\boldsymbol{\Sigma}) \in \Theta_{{\boldsymbol{\vartheta}}_2}$. We assume $\Theta_{{\boldsymbol{\vartheta}}_2}$ is compact.
\begin{assumption} \label{assn_compact}
The parameter space $\Theta_{{\boldsymbol{\vartheta}}_2}$ is compact.
\end{assumption}
Assume $\alpha\in [0,3/4]$ without loss of generality. Then, the null hypothesis $H_0:M=1$ is written as
\[
H_0: \alpha(\boldsymbol{\mu}_1-\boldsymbol{\mu}_2)=0.
\]
For an arbitrary small $\zeta>0$, we partition the parameter space as $\Theta_{{\boldsymbol{\vartheta}}_2} = \Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^1 \cup \Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^2 $, where
\[
\Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^1 =\{{\boldsymbol{\vartheta}}_2\in \Theta_{{\boldsymbol{\vartheta}}_2} : |\boldsymbol{\mu}_{1}- \boldsymbol{\mu}_{2}|\leq \zeta \}\ \ \text{and}\ \ \Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^2 =\{{\boldsymbol{\vartheta}}_2\in \Theta_{{\boldsymbol{\vartheta}}_2}: |\boldsymbol{\mu}_{1}- \boldsymbol{\mu}_{2}|  \geq \zeta \}.
\]
Let $L_n({\boldsymbol{\vartheta}}_2): = \sum_{i = 1}^n \log f_2( \boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_2)$ denote the log-likelihood function, and define the two-component MLE by $\widehat{\boldsymbol{\vartheta}}_2 :=\arg\max_{ {\boldsymbol{\vartheta}}_2 \in \Theta_{{\boldsymbol{\vartheta}}_2}} L_n({\boldsymbol{\vartheta}}_2)$. Define the restricted two-component MLE by $\widehat{\boldsymbol{\vartheta}}_{2,\zeta}^j :=\arg\max_{ {\boldsymbol{\vartheta}}_2^j \in \Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^j } L_n({\boldsymbol{\vartheta}}_2)$ for $j=1,2$ so that $L_n(\widehat {\boldsymbol{\vartheta}}_2) = \max\{ L_n(\widehat {\boldsymbol{\vartheta}}_{2,\zeta}^1 ),L_n(\widehat {\boldsymbol{\vartheta}}_{2,\zeta}^2 ) \}$. Let $(\widehat{\boldsymbol{\gamma}}_0,\widehat{{\mu}}_0,\widehat{{\sigma}}^2_0)$ denote the one-component MLE that maximizes the one-component log-likelihood function $L_{0,n}(\boldsymbol{\gamma},{ {\mu}},\boldsymbol{\sigma}) := \sum_{i=1}^n \log f ( \boldsymbol{X}_i|\boldsymbol{Z}_i;\boldsymbol{\gamma},{{\mu}}, {\sigma^2})$. Define the LRTS for testing $H_0$ as $LR_{n}:= 2\{L_n(\widehat {\boldsymbol{\vartheta}}_2) - L_{0,n}(\widehat{\boldsymbol{\gamma}}_0,\widehat{{\mu}}_0,\widehat{{\sigma}}^2_0)\}=\max \{LR_{n,\zeta}^1 ,LR_{n,\zeta}^2 \}$, where $LR_{n,\zeta}^j := 2\{L_n(\widehat{\boldsymbol{\vartheta}}_{2,\zeta}^j ) - L_{0,n}(\widehat{\boldsymbol{\gamma}}_0,\widehat{{\mu}}_0,\widehat{{\sigma}}^2_0)\}$.

In the following, we derive the asymptotic null distribution of $LR_{n,\zeta}^1$, $LR_{n,\zeta}^2$, and $LR_{n}$ by using a similar approach to the heteroscedastic case. We introduce a reparameterization that extracts the direction in which the Fisher information matrix is singular and approximate the log-likelihood in terms of the polynomials of the reparameterized parameters.

\subsubsection{Asymptotic distribution of $LR_{n,\zeta}^1$}
In this section, we derive the asymptotic distribution of $LR_{n,\zeta}^1$. Let $\boldsymbol{v}=\boldsymbol{w}(\boldsymbol{\Sigma})$ and consider the following one-to-one mapping between $(\boldsymbol{\mu}_1,\boldsymbol{\mu}_2,\boldsymbol{v})$ and the reparameterized parameter $(\boldsymbol{\lambda} ,\boldsymbol{\nu}_{\boldsymbol\mu},\boldsymbol{\nu}_{\boldsymbol{v}})$:
\begin{equation} \label{repara2-homo}
\begin{pmatrix}
{\boldsymbol{\mu}}_1\\
{\boldsymbol{\mu}}_2\\
\boldsymbol{v}
\end{pmatrix}
 =
\begin{pmatrix}
\boldsymbol{\nu}_{\boldsymbol\mu} + (1-\alpha) \boldsymbol{\lambda} \\
\boldsymbol{\nu}_{\boldsymbol\mu} -\alpha \boldsymbol{\lambda} \\
\boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(1-\alpha) \boldsymbol{w}(\boldsymbol{\lambda} \boldsymbol{\lambda} ^{\top})
\end{pmatrix}.
\end{equation}
In the reparameterized model, the density is given by
\begin{equation} \label{loglike-homo}
\begin{aligned}
g(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\psi},\alpha) & = \alpha f_v\left(\boldsymbol{x}\middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu}+(1-\alpha)\boldsymbol{\lambda}, \boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(1-\alpha) \boldsymbol{w}(\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top})\right) \\
& \quad + (1 - \alpha) f_v \left(\boldsymbol{x} \middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu} -\alpha\boldsymbol{\lambda},\boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(1-\alpha) \boldsymbol{w}(\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top}) \right).
\end{aligned}
\end{equation}

We partition $\boldsymbol{\psi}$ as $\boldsymbol{\psi} = (\boldsymbol{\eta}^{\top},\boldsymbol{\lambda}^{\top})^{\top}$, where $\boldsymbol{\eta}: = (\boldsymbol{\gamma}^{\top},\boldsymbol{\nu}_{\boldsymbol\mu}^{\top},\boldsymbol{\nu}_{\boldsymbol{v}}^{\top})^{\top} \in \Theta_{\boldsymbol{\eta}}$  and $\boldsymbol{\lambda} \in \Theta_{\boldsymbol{\lambda}}$. Denote the true values of $\boldsymbol{\eta}$, $\boldsymbol{\lambda}$, and $\boldsymbol{\psi}$ under $H_{0}$ by $\boldsymbol{\eta}^*: = ((\boldsymbol{\gamma}^*)^{\top},({\boldsymbol{\mu}}^*)^{\top},(\boldsymbol{v}^{*})^{\top})^{\top}$, $\boldsymbol{\lambda}^*: = \boldsymbol{0}$, and $\boldsymbol{\psi}^* = ((\boldsymbol{\eta}^*)^{\top}, \boldsymbol{0}^{\top})^{\top}$, respectively. The first derivative of (\ref{loglike-homo}) w.r.t.\ $\boldsymbol{\eta}$ under $\boldsymbol{\psi} = \boldsymbol{\psi}^*$ is given by (\ref{dvareta}) and the first and second derivatives of $g(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\psi},\alpha)$ w.r.t.\ $\boldsymbol{\lambda}$ become zero when evaluated at $\boldsymbol{\psi}=\boldsymbol{\psi}^*$. Consequently, the information on $\boldsymbol{\lambda}$ is provided by the third and fourth derivatives w.r.t.\ $\boldsymbol{\lambda}$.

Define the score vector $\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z})$ as
\begin{equation} \label{score_defn-homo}
\begin{aligned}
\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z}) & :=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\eta}}\\
\boldsymbol{s}_{\boldsymbol{\lambda}}
\end{pmatrix} :=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\eta}}\\
\boldsymbol{s}_{\boldsymbol{\mu^3}}\\
\boldsymbol{s}_{\boldsymbol{\mu}^4}
\end{pmatrix}
\quad \text{with }\
\underset{(d \times 1)}{\boldsymbol{s}_{\boldsymbol{\eta}}} := \frac{\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top},\boldsymbol{v}^{\top})^{\top}}f^*_v }{ f^*_v}, \\
\underset{(d_{\mu^3} \times 1)}{\boldsymbol{s}_{\boldsymbol{\mu}^3}} &:= \left\{\frac{\nabla_{\mu_i \mu_j \mu_k} f^*_v}{3! f^*_v} \right\}_{1 \leq i \leq j \leq k \leq d}, \quad
\underset{(d_{\mu^4} \times 1)}{\boldsymbol{s}_{\boldsymbol{\mu}^4}} := \left\{ \frac{\nabla_{\mu_i \mu_j \mu_k \mu_\ell} f^*_v }{ 4!f^*_v} \right\}_{1 \leq i \leq j \leq k \leq \ell \leq d} ,
\end{aligned}
\end{equation}
where we suppress the dependence of $(\boldsymbol{s}_{\boldsymbol{\eta}}, \boldsymbol{s}_{\boldsymbol{\mu}^3}, \boldsymbol{s}_{\boldsymbol{\mu}^4})$ on $(\boldsymbol{x},\boldsymbol{z})$. Collect the relevant reparameterized parameters as
\begin{equation} \label{tpsi_defn-homo}
\boldsymbol{t}(\boldsymbol{\psi},\alpha) :=
\begin{pmatrix}
\boldsymbol{\eta}-\boldsymbol{\eta}^*\\
\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)
\end{pmatrix}
:=
\begin{pmatrix}
\boldsymbol{\eta}-\boldsymbol{\eta}^*\\
 \alpha(1-\alpha)(1-2\alpha)\boldsymbol{\lambda}_{\boldsymbol{\mu}^3}\\
\alpha(1-\alpha)(1-6\alpha+6\alpha^2)\boldsymbol{\lambda}_{\boldsymbol{\mu}^4}
\end{pmatrix},
\end{equation}
where $\underset{(d_{\mu^3} \times 1)}{\boldsymbol{\lambda}_{\boldsymbol{\mu}^3}} = \{(\boldsymbol{\lambda}_{\boldsymbol{\mu}^3})_{ijk}\}_{1 \leq i \leq j \leq k \leq d}$ with $(\boldsymbol{\lambda}_{\boldsymbol{\mu }^3})_{ijk} := \sum_{(t_1,t_2,t_3) \in p(i,j,k)} \lambda_{t_1} \lambda_{t_2} \lambda_{t_3}$, where $\sum_{(t_1,t_2,t_3) \in p(i,j,k)}$ denotes the sum over all distinct permutations of $(i,j,k)$ to $(t_1,t_2,t_3)$, and $\boldsymbol{\lambda}_{\boldsymbol{\mu}^4}$ is given in (\ref{lambda_muv_defn}).

The third and fourth order derivatives of the density ratio w.r.t. $\boldsymbol{\lambda}$ are given by $\alpha(1-\alpha)(1-2\alpha)\boldsymbol{\lambda}_{\boldsymbol{\mu}^3}^{\top} \boldsymbol{s}_{\boldsymbol{\mu^3}}$ and $\alpha(1-\alpha)(1-6\alpha+6\alpha^2)\boldsymbol{\lambda}_{\boldsymbol{\mu}^4}^{\top}\boldsymbol{s}_{\boldsymbol{\mu}^4}$, respectively. When $\alpha$ is bounded away from $1/2$, the third order derivative identifies $\boldsymbol{\lambda}$ because the third order derivative dominates the fourth order derivative as $\boldsymbol{\lambda}\to 0$.
When $\alpha=1/2$, the third order derivative is identically equal to zero and the fourth order derivative identifies $\boldsymbol{\lambda}$. When $\alpha$ is in the neighborhood of $1/2$ such that $1-2\alpha \propto \boldsymbol{\lambda}$, the third and fourth order derivatives jointly identify $\boldsymbol{\lambda}$.

Accordingly, we characterize  the limit of possible values of $\sqrt{n}\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)$ defined in (\ref{tpsi_defn-homo}) as $n\rightarrow\infty$ by the following two sets:
\begin{equation} \label{Lambda-e-homo}
\begin{aligned}
\Lambda_{\boldsymbol{\lambda} }^1 &:= \left\{\left( \boldsymbol{t}_{\boldsymbol{\mu}^3}^{\top}, \boldsymbol{t}_{\boldsymbol{\mu}^4}^{\top} \right)^{\top}\in \mathbb{R}^{d_{\mu^3}+d_{\mu^4}}:\ \boldsymbol{t}_{\boldsymbol{\mu}^3}= \boldsymbol{\lambda_{\mu^3}},\ \boldsymbol{t}_{\boldsymbol{\mu}^4}= \boldsymbol{0} \text{ for some $ \boldsymbol{\lambda} \in \mathbb{R}^{d}$ }\right\}, \\
\Lambda_{\boldsymbol{\lambda} }^2 &:= \left\{\left( \boldsymbol{t}_{\boldsymbol{\mu}^3}^{\top}, \boldsymbol{t}_{\boldsymbol{\mu}^4}^{\top} \right)^{\top}\in \mathbb{R}^{d_{\mu^3}+d_{\mu^4}}:\ \boldsymbol{t}_{\boldsymbol{\mu}^3}=c \boldsymbol{\lambda_{\mu^3}},\ \boldsymbol{t}_{\boldsymbol{\mu}^4}=- \boldsymbol{\lambda_{\mu^4}}\text{ for some $(\boldsymbol{\lambda}^{\top},c)^{\top} \in \mathbb{R}^{d+1}$ }\right\},
\end{aligned}
\end{equation}
where $\Lambda_{\boldsymbol{\lambda} }^1$ represents the case when $\alpha$ is bounded away from $1/2$ while, by choosing different values of $c$,  $\Lambda_{\boldsymbol{\lambda} }^2$ represents both  cases when $\alpha=1/2$  and when   $\alpha$ is in the neighborhood of $1/2$.

Define $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^j $ by
\begin{equation} \label{t-lambda-homo}
r(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^j ) = \inf_{\boldsymbol{t}_{\boldsymbol{\lambda}} \in \Lambda_{\boldsymbol{\lambda}}^j } r(\boldsymbol{t}_{\boldsymbol{\lambda}}), \quad r(\boldsymbol{t}_{\boldsymbol{\lambda}}) := (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}),
\end{equation}
where $\Lambda_{\boldsymbol{\lambda} }^j$ for $j=1,2$ is defined in (\ref{Lambda-e-homo}) while $\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}$ and $\boldsymbol{Z}_{\boldsymbol{\lambda}}$ are defined by
\begin{align*}
& \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}:=\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}}-\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda\eta}}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}}^{-1}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta\lambda}} \quad\text{and}\quad \boldsymbol{Z}_{\boldsymbol{\lambda}}:=(\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})^{-1} \boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}   \quad\text{with}\\
&\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}} := E[\boldsymbol{s}_{\boldsymbol{\eta}}\boldsymbol{s}_{\boldsymbol{\eta}}^{\top}], \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}} := E[\boldsymbol{s}_{\boldsymbol{\lambda}}\boldsymbol{s}_{\boldsymbol{\lambda}}^{\top}], \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda\eta} } := E[\boldsymbol{s}_{\boldsymbol{\lambda}}\boldsymbol{s}_{\boldsymbol{\eta}}^{\top}],\quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\eta \lambda}} := \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda \eta}}^{\top},
\end{align*}
given $(\boldsymbol{s}_{\boldsymbol{\eta}}, \boldsymbol{s}_{\boldsymbol{\lambda}})$ in (\ref{score_defn-homo}), where $\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \sim N(0,\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})$.

The following proposition establishes the asymptotic null distribution of $LR_{n,\zeta}^1$.
\begin{proposition} \label{P-LR-N1-homo-1}
Suppose that Assumptions \ref{assn_consis} and \ref{A-taylor1} hold and $\boldsymbol{\mathcal{I}} := E[\boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z})\boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z})^{\top}]$ is finite and nonsingular. Then, under the null hypothesis of $H_{01}: \mu_1=\mu_2$, for any $\zeta>0$, $LR_{n,\zeta}^1 \rightarrow_d \max\left\{ (\widehat{\boldsymbol{t}}^1_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}^1_{\boldsymbol{\lambda}}, (\widehat{\boldsymbol{t}}^2_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}^2_{\boldsymbol{\lambda}} \right\}$.
\end{proposition}
\begin{example}
When $d=1$ with $\boldsymbol{\lambda}=\lambda$, we have $\boldsymbol{s_{\mu^3}} = {\nabla_{\mu^3 } f^*_v}/{3! f^*_v}$, $\boldsymbol{s_{\mu^4}} = {\nabla_{\mu^4} f^*_v}/{4! f^*_v}$, and the possible values of $\sqrt{n}\boldsymbol{t}_{\boldsymbol{\lambda}}(\boldsymbol{\lambda},\alpha)=\sqrt{n}\alpha(1-\alpha) \left( (1-2\alpha)\lambda^3, (1-6\alpha+6\alpha^2)\lambda^4\right)^{\top}$ as $n\to \infty$ are given by $\Lambda_{\lambda}^1\cup \Lambda_{\lambda}^2=\mathbb{R}\times \mathbb{R}_-$. In this case, $LR_{n,\zeta}^1 \rightarrow_d (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}} $ with $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}} $ defined by
$r(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}} ) = \inf_{\boldsymbol{t}_{\boldsymbol{\lambda}} \in \mathbb{R}\times \mathbb{R}_-} r(\boldsymbol{t}_{\boldsymbol{\lambda}})$.
 When $d=2$ with $\boldsymbol{\lambda}=(\lambda_1,\lambda_2)^{\top}$, we have   $\Lambda_{\boldsymbol{\lambda} }^1= \left\{(\boldsymbol{\lambda_{\mu^3}}^{\top},\boldsymbol{0}^{\top})^{\top}: (\lambda_1,\lambda_2)^{\top} \in \mathbb{R}^2\right\}$ and $
\Lambda_{\boldsymbol{\lambda} }^2 := \left\{(c\boldsymbol{\lambda_{\mu^3}}^{\top},-\boldsymbol{\lambda_{\mu^4}}^{\top})^{\top}: (\lambda_1,\lambda_2,c)^{\top} \in \mathbb{R}^3\right\}$ with $\boldsymbol{\lambda_{\mu^3}} = (\lambda_1^3,3\lambda_1^2\lambda_2,3\lambda_1\lambda_2^2,\lambda_2^3)^{\top}$ and $\boldsymbol{\lambda_{\mu^4}} = (\lambda_1^4,4\lambda_1^3\lambda_2,6\lambda_1^2\lambda_2^2,4\lambda_1\lambda_2^3,\lambda_2^4)^{\top}$.
\end{example}




\subsubsection{Asymptotic distribution of $LR_{n,\zeta}^2$}

This section derives the asymptotic distribution of $LR_{n,\zeta}^2$. We use the reparameterization (\ref{repara2-homo}) but collect the reparameterized parameters into $\boldsymbol{\phi}:=(\boldsymbol{\eta}^{\top},\alpha)^{\top}$ and $\boldsymbol{\lambda}$, where $\boldsymbol{\eta}: = (\boldsymbol{\gamma}^{\top},\boldsymbol{\nu}_{\boldsymbol\mu}^{\top},\boldsymbol{\nu}_{\boldsymbol{v}}^{\top})^{\top}$. Let the resulting density be
\begin{equation} \label{loglike-homo-2}
\begin{aligned}
h(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\phi},\boldsymbol{\lambda}) & := \alpha f_v\left(\boldsymbol{x}\middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu}+(1-\alpha)\boldsymbol{\lambda}, \boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(1-\alpha) \boldsymbol{w}(\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top})\right) \\
& \quad + (1 - \alpha) f_v \left(\boldsymbol{x} \middle|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\nu}_{\boldsymbol\mu} -\alpha\boldsymbol{\lambda},\boldsymbol{\nu}_{\boldsymbol{v}} - \alpha(1-\alpha) \boldsymbol{w}(\boldsymbol{\lambda}\boldsymbol{\lambda}^{\top}) \right).
\end{aligned}
\end{equation}
The right hand side of (\ref{loglike-homo-2}) is the same as that of (\ref{loglike-homo}). When we restrict the parameter space to $\Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^2$, the reparameterized density $h(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\phi},\boldsymbol{\lambda})$ in (\ref{loglike-homo-2}) becomes the one-component density if and only if $\alpha =0$. Furthermore, $\boldsymbol{\lambda}$ is not identified when $\alpha=0$. Denote the true value of $\boldsymbol{\phi}$ under $H_{0}$ by $\boldsymbol{\phi}^* = ((\boldsymbol{\eta}^*)^{\top}, 0)^{\top}$, where $\boldsymbol{\eta}^* := ((\boldsymbol{\gamma}^*)^{\top},({\boldsymbol{\mu}}^*)^{\top},(\boldsymbol{v}^{*})^{\top})^{\top}$. The MLE of $\boldsymbol{\phi}$ under the restriction $\boldsymbol{\vartheta}_2\in\Theta_{{\boldsymbol{\vartheta}}_2,\zeta}^2$ converges to $\boldsymbol{\phi}^*$ in probability.

Define $f_v^*(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\lambda}):=f_v\left(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*, \boldsymbol\mu^*+ \boldsymbol{\lambda}, \boldsymbol{v}^*\right)$ so that $f^*_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{0}) = f^*_v$ and $\nabla f^*_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{0}) = \nabla f^*_v$. Define the score vectors $\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z};\boldsymbol{\lambda})$ indexed by $\boldsymbol{\lambda}$ as
\begin{equation} \label{score_defn-homo-2}
\boldsymbol{s}(\boldsymbol{x},\boldsymbol{z};\boldsymbol{\lambda}) :=
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\eta}}\\
{s}_{\alpha}(\boldsymbol{\lambda})
\end{pmatrix},
\end{equation}
where $\boldsymbol{s}_{\boldsymbol{\eta}} =\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top},\boldsymbol{v}^{\top})^{\top}}f_v^*/f_v^*$ as defined in (\ref{score_defn-homo}) and
\begin{equation} \label{score_defn-homo-3}
{s}_{\alpha}(\boldsymbol{\lambda}) :=
\frac{f_v^*(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\lambda}) -f_v^* - \nabla_{\boldsymbol{\mu}^{\top}} f_v^* \boldsymbol{\lambda} - \nabla_{\boldsymbol{v}^{\top}} f_v^* \boldsymbol{\lambda}_{\boldsymbol{\mu}^2} }{ |\boldsymbol{\lambda}|^3 f_v^*},
\end{equation}
where $\underset{(d_{\mu^2} \times 1)}{\boldsymbol{\lambda}_{\boldsymbol{\mu}^2}} := \{(\boldsymbol{\lambda}_{\boldsymbol{\mu}^2})_{ij}\}_{1 \leq i \leq j \leq d}$ with $(\boldsymbol{\lambda}_{\boldsymbol{\mu }^2})_{ij}:=\lambda_{ii}^2$ if $i=j$ and $2\lambda_{ij}$ if $i \neq j$. The division by $|\boldsymbol{\lambda}|^3$ is necessary to define $s_\alpha(\boldsymbol{\lambda})$ here because, if we were to define $s_\alpha(\boldsymbol{\lambda})$ as $(f_v^*(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\lambda}) -f_v^* - \nabla_{\boldsymbol{\mu}^{\top}} f_v^* \boldsymbol{\lambda} - \nabla_{\boldsymbol{v}^{\top}} f_v^* \boldsymbol{\lambda}_{\boldsymbol{\mu}^2})/ f_v^*$, then we have $s_\alpha(\boldsymbol{\lambda})\to 0$ as $\boldsymbol{\lambda}\to 0$, invalidating the approximation when $\boldsymbol{\lambda}$ is close to zero.

Collect the relevant reparameterized parameters as
\begin{equation} \label{tpsi_defn-homo-2}
\boldsymbol{t} (\boldsymbol{\phi},\boldsymbol\lambda)
:=
\begin{pmatrix}
\boldsymbol{t}_{\boldsymbol{\eta}}(\boldsymbol\lambda)\\
t_\alpha(\boldsymbol{\lambda})
\end{pmatrix}\quad\text{with}\quad
\boldsymbol{t}_{\boldsymbol{\eta}}(\boldsymbol\lambda):=
\begin{pmatrix}
\boldsymbol{\gamma}-\boldsymbol{\gamma}^*\\
\boldsymbol{\nu_\mu}-\boldsymbol{\mu}^* \\
\boldsymbol{\nu_v}-\boldsymbol{v}^*
\end{pmatrix}\quad \text{and} \quad
t_\alpha(\boldsymbol{\lambda}):=\alpha |\boldsymbol{\lambda}|^3.
\end{equation}

In (\ref{score_defn-homo-3}), $s_\alpha(\boldsymbol{\lambda})$ is non-degenerate and not perfectly correlated with $\boldsymbol{s}_{\boldsymbol{\eta}}$ even when $\boldsymbol{\lambda} \to \boldsymbol{0}$. With $(\boldsymbol{s}_{\boldsymbol{\eta}}, {s}_{\alpha}(\boldsymbol{\lambda}))$ defined in (\ref{score_defn-homo-2}), define
\begin{equation} \label{I_lambda-homo-2}
\begin{aligned}
&\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}} := E[\boldsymbol{s}_{\boldsymbol{\eta}}(\boldsymbol{s}_{\boldsymbol{\eta}})^{\top}], \quad \boldsymbol{\mathcal{I}}_{\alpha\boldsymbol{\eta} }(\boldsymbol{\lambda}) := E[{s}_{\alpha}(\boldsymbol{\lambda}) \boldsymbol{s}_{\boldsymbol{\eta}} ^{\top}], \quad \boldsymbol{\mathcal{I}}_{\boldsymbol{\eta }\alpha}(\boldsymbol{\lambda}) := (\boldsymbol{\mathcal{I}}_{\alpha \boldsymbol{\eta}}(\boldsymbol{\lambda}))^{\top},\\
& {\mathcal{I}}_{\alpha}(\boldsymbol{\lambda}_1,\boldsymbol{\lambda}_2) := E[{s}_{\alpha}(\boldsymbol{\lambda}_1) {s}_{\alpha}(\boldsymbol{\lambda}_2) ],\quad {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}_1,\boldsymbol{\lambda}_2):= {\mathcal{I}}_{\alpha}(\boldsymbol{\lambda}_1,\boldsymbol{\lambda}_2)-\boldsymbol{\mathcal{I}}_{\alpha\boldsymbol{\eta}}(\boldsymbol{\lambda}_1)(\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}})^{-1}\boldsymbol{\mathcal{I}}_{\boldsymbol{\eta}\alpha}(\boldsymbol{\lambda}_2), \\
& {Z}_{\alpha}(\boldsymbol{\lambda}):=( {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}, \boldsymbol{\lambda}))^{-1} {G}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}),
\end{aligned}
\end{equation}
where ${G}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda})$ is a mean zero Gaussian process indexed by $\boldsymbol{\lambda}$ with $\text{Cov}({G}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}_1),{G}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}_2)) = {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}_1,\boldsymbol{\lambda}_2)$.
Define $\widehat{t}_\alpha(\boldsymbol{\lambda})$ by
\begin{align} \label{t-lambda-homo-2}
r(\widehat{t}_\alpha(\boldsymbol{\lambda}))= \inf_{ t_\alpha \geq 0 }r( t_\alpha),\quad
 r( t_\alpha) := (t_\alpha - {Z}_{{\alpha}}(\boldsymbol{\lambda}))^2 {\mathcal{I}} _{{\alpha}.\boldsymbol{\eta}}(\boldsymbol{\lambda}, \boldsymbol{\lambda}).
 \end{align}
The following proposition establishes the asymptotic null distribution of $LR_{n,\zeta}^{2}$.  Define $\boldsymbol{\mathcal{I}}(\boldsymbol{\lambda}) := E[\boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z};\boldsymbol{\lambda}) \boldsymbol{s}(\boldsymbol{X},\boldsymbol{Z};\boldsymbol{\lambda})^{\top} ]$.
\begin{proposition} \label{P-LR-N1-homo-2}
Suppose that Assumptions \ref{assn_consis}, \ref{A-taylor1}, and \ref{assn_compact} hold and $0 < \inf_{\Theta_{\boldsymbol{\lambda}} \setminus \{\boldsymbol{0}\}} \lambda_{\min}(\boldsymbol{\mathcal{I}}(\boldsymbol{\lambda})) \leq \sup_{\Theta_{\boldsymbol{\lambda}} \setminus \{\boldsymbol{0}\}} \lambda_{\max}(\boldsymbol{\mathcal{I}}(\boldsymbol{\lambda})) < \infty$. Then, under the null hypothesis of $H_{0}: M=1$, for any $\zeta>0$, $LR_{n,\zeta}^{2} \rightarrow_d \sup_{\Theta_{\boldsymbol{\lambda}} \cap \{|\boldsymbol{\lambda}|\geq \zeta\} }\ (\widehat{t}_\alpha(\boldsymbol{\lambda}) )^2 {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}, \boldsymbol{\lambda})$.
\end{proposition}

\subsubsection{Testing $H_{0}: M=1$}

The following proposition derives the asymptotic distribution of $LR_{n}$.
The proof is omitted because it is a straightforward consequence of Propositions \ref{P-LR-N1-homo-1} and \ref{P-LR-N1-homo-2} in view of $LR_{n}=\lim_{\zeta\to 0}\max \{LR_{n,\zeta}^1 ,LR_{n,\zeta}^2 \}$.

\begin{proposition}\label{P-LR-N1-homo}
Suppose that Assumptions of Propositions \ref{P-LR-N1-homo-1} and \ref{P-LR-N1-homo-2} hold. Then, under the null hypothesis of $H_{0}: M=1$,
\[
LR_{n} \rightarrow_d \max\left\{ (\widehat{\boldsymbol{t}}^1_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}^1_{\boldsymbol{\lambda}}, (\widehat{\boldsymbol{t}}^2_{\boldsymbol{\lambda}})^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}} \widehat{\boldsymbol{t}}^2_{\boldsymbol{\lambda}}, \ \sup_{\Theta_{\boldsymbol{\lambda}} \setminus \{\boldsymbol{0}\} }\ (\widehat{t}_\alpha(\boldsymbol{\lambda}) )^2 {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}, \boldsymbol{\lambda}) \right\},
\]
where $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^j $ for $j=1,2$ is defined in (\ref{t-lambda-homo}) and $\widehat{t}_\alpha(\boldsymbol{\lambda})$ is defined in (\ref{t-lambda-homo-2}).
\end{proposition}
This result generalizes Theorem 2 of \citet{chenchen03sinica}, who derive the asymptotic distribution of the LRTS in the univariate case.


\subsection{Likelihood ratio test of $H_0: M = M_0$ against $H_A: M = M_0 + 1$ for $M_0\geq 2$}

We consider a random sample $\{\boldsymbol{X}_i,\boldsymbol{Z}_i\}_{i = 1}^n$ generated from the following $M_0$-component $d$-variate normal mixture density model with common variance:
\begin{equation}
f_{M_0}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0}^*):=\sum_{j=1}^{M_0} \alpha_{j}^{*} f(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\gamma}}^*, \boldsymbol{\mu}_{j}^{*},\boldsymbol{\Sigma}^{*}), \label{true_model-homo}
\end{equation}
where ${\boldsymbol{\vartheta}}_{M_0}^*=(\alpha_{1}^*,\ldots,\alpha_{M_0 - 1}^*,{\boldsymbol{\gamma}}^*,\boldsymbol{\mu}_1^*,\ldots,\boldsymbol{\mu}_{M}^*,\boldsymbol{\Sigma}^{*} )$ and $\alpha_j^*>0$. We assume $\boldsymbol{\mu}_{1}^* <\ldots< \boldsymbol{\mu}_{M_0}^*$ for identification. The corresponding density of an $(M_0+1)$-component mixture model is given by
\begin{equation}
f_{M_0+1}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0+1}):=\sum_{j=1}^{M_0+1}\alpha_j f(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma},\boldsymbol{\mu}_j,\boldsymbol{\Sigma} ),\label{fitted_model-homo}
\end{equation}
where ${\boldsymbol{\vartheta}}_{M_0+1} = (\alpha_1,\ldots,\alpha_{M_0},\boldsymbol{\gamma},\boldsymbol{\mu}_1,\ldots.,\boldsymbol{\mu}_{M_0+1},\boldsymbol{\Sigma})$. Partition the null hypothesis as $H_0 = \cup_{m=1}^{M_0} H_{0,m}$ with $H_{0,m}: \alpha_m(\boldsymbol{\mu}_{m} - \boldsymbol{\mu}_{m + 1}) = 0$.

Define the LRTS for testing $H_{01}$ as
\[
 LR_{n}^{M_0} := \max_{{\boldsymbol{\vartheta}}_{M_0+1}\in \Theta_{{\boldsymbol{\vartheta}}_{M_0+1}}}2\{ L_n({\boldsymbol{\vartheta}}_{M_0+1})- L_{0,n}(\widehat{{\boldsymbol{\vartheta}}}_{M_0})\},
\]
where $L_n({\boldsymbol{\vartheta}}_{M_0 + 1}):=\sum_{i = 1}^n \log f_{M_0 + 1}(\boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_{M_0 + 1})$, $L_{0,n}({\boldsymbol{\vartheta}}_{M_0})=\sum_{i = 1}^n \log f_{M_0}(\boldsymbol{X}_i|\boldsymbol{Z}_i;{\boldsymbol{\vartheta}}_{M_0})$, and $\widehat{\boldsymbol{\vartheta}}_{M_0}=\mathop{\arg \max}_{{\boldsymbol{\vartheta}}_{M_0}\in\Theta_{{\boldsymbol{\vartheta}}_{M_0}}}L_{0,n}({\boldsymbol{\vartheta}}_{M_0})$ for the densities (\ref{true_model-homo})--(\ref{fitted_model-homo}). Define $\widetilde{\boldsymbol{\lambda}}:=((\boldsymbol{\lambda}_1)^{\top},\ldots,(\boldsymbol{\lambda}_{M_0})^{\top})^{\top} \in \Theta_{\widetilde{\boldsymbol{\lambda}}}$ with $\boldsymbol{\lambda}_m \in \Theta_{\boldsymbol{\lambda}_m}$.  Collect the score vector for testing $H_{0,1},\ldots,H_{0,M_0}$   as
\begin{equation}
\begin{aligned}
&\widetilde{\boldsymbol{s}}(\boldsymbol{x},\boldsymbol{z}) :=
\begin{pmatrix}
\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}} \\ \widetilde{\boldsymbol{s}}_{\boldsymbol{\lambda}}
\end{pmatrix}\quad \text{and}\quad
\bar{\boldsymbol{s}}(\boldsymbol{x},\boldsymbol{z};\widetilde{\boldsymbol{\lambda}}):=
\begin{pmatrix}
\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}} \\ \bar{\boldsymbol{s}}_{\boldsymbol\alpha}(\tilde{\boldsymbol{\lambda}})
\end{pmatrix},
\ \text{ where }\\
&\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}}:= \left(\begin{array}{c} {\boldsymbol{s}}_{\boldsymbol{\alpha}} \\ {\boldsymbol{s}}_{{ (\boldsymbol{\gamma},\boldsymbol{\mu}, \boldsymbol{v}) }}
 \end{array}\right),\quad \widetilde{\boldsymbol{s}}_{\boldsymbol{\lambda}}: =
\begin{pmatrix}
\boldsymbol{s}_{\boldsymbol{\mu}^3}^1 \\
\boldsymbol{s}_{\boldsymbol{\mu}^4}^1 \\
 \vdots \\
\boldsymbol{s}_{\boldsymbol{\mu}^3}^{M_0} \\
\boldsymbol{s}_{\boldsymbol{\mu}^4}^{M_0}
\end{pmatrix}, \quad \text{and}\quad
\bar{\boldsymbol{s}}_{\boldsymbol\alpha}(\widetilde{\boldsymbol{\lambda}}): =
\begin{pmatrix}
{s}_{\boldsymbol{\alpha}}^1(\boldsymbol{\lambda}_1)\\
 \vdots \\
{s}_{\boldsymbol{\alpha}}^{M_0}(\boldsymbol{\lambda}_{M_0})
\end{pmatrix},
\end{aligned}
\label{stilde-homo}
\end{equation}
where $\boldsymbol{s}_{\boldsymbol{\mu^3}}^m := \left\{\alpha_{m}^* \nabla_{\mu_i \mu_j \mu_k} f^*_v(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*,\boldsymbol{\mu}_m^{*},\boldsymbol{v}^*) / (3!f_0^*) \right\}_{1 \leq i \leq j \leq k \leq d}$; $\boldsymbol{s}_{\boldsymbol{\alpha}}$, $\boldsymbol{s}_{(\boldsymbol{\gamma},\boldsymbol{\mu}, \boldsymbol{v})}$, and $\boldsymbol{s}_{\boldsymbol{\mu}^4}^m$ are defined similarly to those in (\ref{sh}) but using the density (\ref{true_model-homo}) in place of (\ref{true_model}) with the common value of $v^*$ across components; ${s}_{\alpha}^m(\boldsymbol{\lambda}_m)$ is defined as
\begin{equation*}
{s}_{\alpha}^m(\boldsymbol{\lambda}_m) := \alpha_m^*
\frac{f_v^{m*}(\boldsymbol{\lambda}_m) -f_v^{m*}-\nabla_{\boldsymbol{\mu}} f^{m*}_v \boldsymbol{\lambda}_m - \nabla_{\boldsymbol{v}} f^{m*}_v \boldsymbol{\lambda}_{\boldsymbol{\mu}^2,m} }{ |\boldsymbol{\lambda}_m|^3f^{m*}_v},
\end{equation*}
where $ f_v^{m*}(\boldsymbol{\lambda}_m):=f_v\left(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*, \boldsymbol\mu_m^*+ \boldsymbol{\lambda}_m, \boldsymbol{v}^*\right)$ and $f_v^{m*} :=f_v\left(\boldsymbol{x}|\boldsymbol{z};\boldsymbol{\gamma}^*, \boldsymbol\mu_m^*, \boldsymbol{v}^*\right)$, and $\boldsymbol{\lambda}_{\boldsymbol{\mu}^2,m} $ is defined similarly to $\boldsymbol{\lambda}_{\boldsymbol{\mu}^2}$ but with $\boldsymbol{\lambda}_m$ in place of $\boldsymbol{\lambda}$.


Define $\widetilde{\boldsymbol{\mathcal{I}}}$, $\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta}}$, $\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}\boldsymbol{\eta}}$, $\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\eta\lambda}}$, $\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}}$, $\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}$ similarly to those in (\ref{Itilde}) but using $\widetilde{\boldsymbol{s}}(\boldsymbol{x},\boldsymbol{z})$ defined in (\ref{stilde-homo}) in place of (\ref{stilde}). Let $\widetilde{\boldsymbol{G}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}=((\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^{1})^{\top},\ldots,(\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^{M_0})^\top)^\top \sim N(0,\widetilde{\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}})$, and define ${\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m:=E[\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m (\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m)^\top]$ and $\boldsymbol{Z}_{\boldsymbol{\lambda}}^m:=({\boldsymbol{\mathcal{I}}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m)^{-1}\boldsymbol{G}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m$. For $j=1,2$, define $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^j $ by
\begin{equation*}
r^m(\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^j ) = \inf_{\boldsymbol{t}_{\boldsymbol{\lambda}} \in \Lambda_{\boldsymbol{\lambda} }^j}r^m(\boldsymbol{t}_{\boldsymbol{\lambda}}), \quad r^m(\boldsymbol{t}_{\boldsymbol{\lambda}}) := (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}^m)^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m (\boldsymbol{t}_{\boldsymbol{\lambda}} -\boldsymbol{Z}_{\boldsymbol{\lambda}}^m),
\end{equation*}
where $\Lambda_{\boldsymbol{\lambda}}^j$ is given by (\ref{Lambda-e-homo}).

Define $\widehat{t}_{\alpha,m}(\boldsymbol{\lambda}_m)$ by
\begin{align*}
r^m(\widehat{t}_{\alpha,m} (\boldsymbol{\lambda}_m))= \inf_{ t_\alpha \geq 0 }r^m( t_\alpha),\quad
 r^m( t_\alpha) := (t_\alpha - {Z}_{\boldsymbol{\alpha}}^m(\boldsymbol{\lambda}_m))^2 {\mathcal{I}} _{\boldsymbol{\alpha}.\boldsymbol{\eta}}^m(\boldsymbol{\lambda}_m, \boldsymbol{\lambda}_m),
\end{align*}
where ${\mathcal{I}}_{\boldsymbol{\alpha}.\boldsymbol{\eta}}^m(\boldsymbol{\lambda}_m, \boldsymbol{\lambda}_m)$ and ${Z}_{\boldsymbol{\alpha}}^m(\boldsymbol{\lambda}_m)$ are defined similarly to $ {\mathcal{I}}_{\alpha.\boldsymbol{\eta}}(\boldsymbol{\lambda}, \boldsymbol{\lambda})$ and $Z_{{\alpha}}(\boldsymbol{\lambda})$ in (\ref{I_lambda-homo-2}), respectively, but using $\widetilde{\boldsymbol{s}}_{\boldsymbol{\eta}}$ and $s_{\boldsymbol{\alpha}}^m(\boldsymbol{\lambda}_m)$ in place of ${\boldsymbol{s}}_{\boldsymbol{\eta}}$ and $s_{\alpha}(\boldsymbol{\lambda})$.


\begin{assumption}\label{A-vec-2-homo}  (a) The parameter spaces $\Theta_{\boldsymbol{\vartheta}_{M_0}}$ and $\Theta_{\boldsymbol{\vartheta}_{M_0+1}}$ are compact. (b)  $\widetilde{\boldsymbol{\mathcal{I}}}=E[\widetilde{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z})\widetilde{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z})^{\top}]$ is finite and nonsingular and $0 < \inf_{\Theta_{\widetilde{\boldsymbol{\lambda}}}\setminus \{\boldsymbol{0}\}  } \lambda_{\min}(\bar{\boldsymbol{\mathcal{I}}}(\widetilde{\boldsymbol{\lambda}})) \leq \sup_{\Theta_{\widetilde{\boldsymbol{\lambda}}}\setminus \{\boldsymbol{0}\}} \lambda_{\max}(\bar{\boldsymbol{\mathcal{I}}}(\widetilde{\boldsymbol{\lambda}})) < \infty$, where $\bar{\boldsymbol{\mathcal{I}}}(\widetilde{\boldsymbol{\lambda}}):= E[\bar{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z};\widetilde{\boldsymbol{\lambda}}) (\bar{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z};\widetilde{\boldsymbol{\lambda}}))^{\top}]$ and $\widetilde{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z})$ and $\bar{\boldsymbol{s}}(\boldsymbol{X},\boldsymbol{Z};\widetilde{\boldsymbol{\lambda}})$ are defined in (\ref{stilde-homo}).
\end{assumption}

\begin{proposition} \label{local_lr-2-homo}
Suppose that Assumptions \ref{assn_consis}, \ref{A-taylor1}, and \ref{A-vec-2-homo} hold. Then, under the null hypothesis $H_0: m=M_0$, $LR_{n}^{M_0} \rightarrow_d \max\{v_1,\ldots, v_{M_0}\}$, where
\[
v_m := \max\left\{ (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^1 )^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^1 , \ (\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^2 )^{\top} \boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}^m \widehat{\boldsymbol{t}}_{\boldsymbol{\lambda},m}^2 , \ \sup_{\Theta_{\boldsymbol{\lambda}_m} \setminus \{\boldsymbol{0}\}}\ (\widehat{t}_{\alpha,m}(\boldsymbol{\lambda}_m) )^2 {\mathcal{I}}_{\boldsymbol{\alpha}.\boldsymbol{\eta}}^m(\boldsymbol{\lambda}_m, \boldsymbol{\lambda}_m) \right\}.
\]
\end{proposition}

\section{Simulation} \label{section:simulation }

\subsection{Choice of penalty function} \label{section:penalty}

To apply our EM test, we need to specify the set $\mathcal{T}$, number of iterations $K$, and penalty functions for  $p(\tau)$. Based on our experience, we recommend $\mathcal{T} = \{0.1, 0.3, 0.5\}$ and $K = \{1,2,3\}$.
We set $p(\tau)= \log(2\min\{\tau,1-\tau\})$ as suggested by \citet{chenli09as}. When estimating the model under the null hypothesis and computing $L_{0,n}(\widehat{\boldsymbol{\vartheta}}_{M_0})$, we use the penalty function (\ref{pen_pmle}) and set $a_n = n^{-1/2}$ as recommended by \citet{chentan09jmva}. For the alternative model,  we consider
$a_n = n^{-1/2}$  and $1$ to examine the sensitivity of the rejection frequencies to the choice of $a_n$.

\subsection{Simulation results} \label{section:results}

We examine the type I error rates and powers of the EM test by small simulations using mixtures of bivariate normal distributions. Computation was done using R \citep{R}. The critical values are computed by bootstrap with $399$ and $199$ bootstrap replications when testing $H_0:M=1$ and $H_0:M=2$, respectively. We use $2000$ replications, and the sample sizes are set to $200$ and $400$.

Table \ref{table1} reports the type I error rates of the EM test of $H_0:M = 1$ against the alternative $H_1:M = 2$  under the null hypothesis using two models given at the bottom of Table \ref{table1}. In both models, the EM test statistics give accurate type I errors for $n=200$ and $400$ across two values of $a_n=n^{-1/2}$ and $1$. Table \ref{table3} reports the powers of the EM test  when  $a_n=1$ under three alternative models given in Table \ref{table2}. Comparing the rejection frequency of Model 1 with that of Model 2 or Model 3, the EM test shows higher power as the distance between two component distributions in the alternative model increases in terms of means (Model 2) or variance (Model 3).


Table \ref{table4} reports the type I error rates of the EM test of $H_0:M = 2$ against the alternative $H_1:M = 3$  under the  two null models given at the bottom of   Table \ref{table4}. The EM test gives accurate type I errors across two models, sample sizes, and the values of $a_n$. Table \ref{table6} reports the powers of the EM test of  $H_0:M = 2$ against the alternative $H_1:M = 3$ under the alternative model given in Table \ref{table5}. Overall, the EM test shows good power under finite sample size.


\section{Empirical applications}
The sequential hypothesis testing based on our  EM test provides a useful alternative to the AIC or the BIC in determining the number of components in empirical applications.\footnote{For penalty function in our empirical applications, we set $a_n=n^{-1/2}$ for the null model and $a_n=1$ for the alternative model. }

\subsection{The flea beetles}

The flea beetles data available in R package \textbf{tourr} contains a sample of 74 flea beetles from three species, ``Concinna," ``Heikertingeri," and ``Heptapotamica'' with  21, 31, and 22 observations, respectively.\footnote{The data is originally from \citet{lubischew62bio}.}  Figure \ref{figure1} provides a scatter plot of two physical measurements ``tars1'' and ``aede1,'' which measure the
width of the first joint of the first tarsus in microns (the sum of measurements for both tarsi) and the maximal width of the aedeagus in the fore-part in microns, respectively, for each of three species. We sequentially test the number of components in this data set without utilizing the information on which species each observation is from.
As shown in Table \ref{table7}, the $p$-values of EM test for testing $H_0: M=1$ and $H_0: M=2$ are $0.00$ and $0.01$, suggesting that the number of components is larger than two. On the other hand, the $p$-values of the EM test for testing $H_0: M=3$ are between 0.32 and 0.36; consistent with the actual number of species in this data set, we fail to reject $H_0: M=3$. In contrast, both the AIC and the BIC incorrectly indicate that there is only one component. Table \ref{table8} compares the estimated three-component bivariate normal mixture model in the first panel  with  the single component models estimated from a subsample of each of three species in the second panel, showing that each of estimated three component distributions accurately captures the corresponding species.

\subsection{Analysis of differential gene expression}

A multivariate normal mixture model can be used to find differentially expressed genes by means of the posterior probability that an individual gene is non-differentially expressed. We analyze the rat dataset of 1,176 genes
in middle-ear mucosa of six rat samples, the first two without  pneumococcal middle-ear infection and the latter four with the disease \citep{pan02bio,he06csda}. As in \citet{pan02bio}, the data were normalized by  log-transformation and median centering. Denote the resulting expression levels of gene $i$  of sample $j$ by $x_{ij}$.
We apply finite bivariate normal mixtures to model the sample average expression levels for gene $i$ under the two conditions,  $(z_{i0},z_{i1})=(\sum_{j=1}^2 x_{ij}/2,\sum_{j=3}^6 x_{ij}/4)$ for $i=1,\ldots,1176$.

As shown in Table \ref{table9},  the sequential hypothesis testing based on EM test and the AIC indicate that there are six components; on the other hand, consistent with the result in \cite{he06csda}, the BIC chooses the five-component model.  Table \ref{table10} presents the estimates from the six component model. We classify each pair of gene expression levels into six clusters using their posterior probabilities and plot them in Figure \ref{figure2}.

32 genes classified into cluster 5 show some evidence for differential expression with a mean difference of 0.23. Similarly,  13 genes in cluster 6 demonstrate strong  evidence for differential expression, albeit with large variability. In contrast, the genes in clusters  1--4 show a flat expression pattern, where the observations in each cluster center around the 45 degree line.


\newpage