Testing the Order of Multivariate Normal Mixture Models

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]

{ \singlespacing }

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

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, lindsay95book, mclachlanpeel00book, and 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 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 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 chen95as. As a result, neither Assumption (P1) of dacunha99as nor Assumption 7 of azais09esaim holds, and Assumption 3 of zhuzhang04jrssb is violated, while Corollary 4.1 of liushao03as does not hold in normal mixtures. Heteroscedastic normal mixture models have an additional problem called the infinite Fisher information problem 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 lcm09bm and 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, chenli09as develop an EM test for $M_0 = 1$ against $M_0=2$, and 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 chenli09as to the multivariate case. 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, chenchen03sinica derive the asymptotic distribution of the LRTS. Our results generalize the results in 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, 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.

Heteroscedastic multivariate finite normal mixture models

Denote the density of a $d$-variate normal distribution with mean $\boldsymbol{\mu}+ \boldsymbol{\gamma}^\top\boldsymbol{z}$ and variance $\boldsymbol{\Sigma}$ by

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:

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)). Our objective is to test \[ H_0:\ M = M_0\quad \text{against}\quad H_A: M = M_0 + 1. \]

Likelihood ratio test of $H_0: M = 1$ against $H_A: M = 2$

As a precursor of the EM test developed in Section (ref), 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

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 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 hartigan85book. Therefore, we consider a maximum penalized likelihood estimator (PMLE) introduced by chentan09jmva. Similar to chentan09jmva, we use the following penalty function with $M=2$:

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)$.

Model ((ref)) 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$.

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

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 rotnitzky00bernoulli and kasaharashimotsu15jasa. Collect the unique elements in $\boldsymbol{\Sigma}$ into a $d(d+1)/2$-vector

Define the density of $N(\boldsymbol{\mu},\boldsymbol{\Sigma})$ parameterized in terms of $\boldsymbol{\mu}$ and $\boldsymbol{v}$ as

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

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}})$:

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

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

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)) 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:

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

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

with $b(\alpha): = -(2/3) (\alpha^2 - \alpha + 1)<0$ and

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)), $\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)$,

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)), define

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)) as $n\rightarrow\infty$. Define

For $j=1,2$, define $\widehat{\boldsymbol{t}}_{\boldsymbol{\lambda}}^{j}$ by

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))--((ref)).

The following proposition establishes the asymptotic null distribution of the LRTS.

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}$.

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

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

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

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$ kasaharashimotsu15jasa. Define the set of values of ${\boldsymbol{\vartheta}}_{M_0+1}$ that yields the true density ((ref)) as

Under $H_{0,1m}$, the $(M_0 + 1)$-component model ((ref)) generates the true $M_0$-component density ((ref)) 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

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

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))--((ref)) and the penalty function in ((ref)). We consider the LRTS for testing $H_{01}$ given by

Collect the score vector for testing $H_{0,11},\ldots,H_{0,1M_0}$ into one vector as

where, with $f_0^*:=f_{M_0}(\boldsymbol{x}|\boldsymbol{z};{\boldsymbol{\vartheta}}_{M_0}^*)$ and for $m=1,\ldots,M_0$,

Define

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

where $\Lambda_{\boldsymbol{\lambda} }^{j}$ is defined in ((ref)). 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.

EM test

Implementing the likelihood ratio test in Section (ref) 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)) 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

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 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

In an M-step, update $\tau$ and $\boldsymbol{\alpha}$ by

and update $\boldsymbol{\mu}_j$ and $\boldsymbol{\Sigma}_j$ for $j=1,\ldots,M_0+1$ by

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 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

where $\widehat{\boldsymbol{\vartheta}}_{M_0}$ and $L_{0,n}({\boldsymbol{\vartheta}}_{M_0})$ are defined in ((ref)).

Finally, with a pre-specified number $K$, define the 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 EM test statistic is defined as the maximum of $M_0$ local EM test statistics:

The following proposition shows that for any finite $K$, the EM test statistic is asymptotically equivalent to the penalized LRTS for testing $H_{01}$.

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)),

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)), 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.

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}$.

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.

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.

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.

Likelihood ratio test of $H_0: M = 1$ against $H_A: M = 2$

Consider a two-component normal mixture density function with common variance:

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.

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.

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}})$:

In the reparameterized model, the density is given by

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)) w.r.t.\ $\boldsymbol{\eta}$ under $\boldsymbol{\psi} = \boldsymbol{\psi}^*$ is given by ((ref)) 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

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

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)).

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)) as $n\rightarrow\infty$ by the following two sets:

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

where $\Lambda_{\boldsymbol{\lambda} }^j$ for $j=1,2$ is defined in ((ref)) while $\boldsymbol{\mathcal{I}}_{\boldsymbol{\lambda}.\boldsymbol{\eta}}$ and $\boldsymbol{Z}_{\boldsymbol{\lambda}}$ are defined by

given $(\boldsymbol{s}_{\boldsymbol{\eta}}, \boldsymbol{s}_{\boldsymbol{\lambda}})$ in ((ref)), 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$.

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

This section derives the asymptotic distribution of $LR_{n,\zeta}^2$. We use the reparameterization ((ref)) 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

The right hand side of ((ref)) is the same as that of ((ref)). 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)) 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

where $\boldsymbol{s}_{\boldsymbol{\eta}} =\nabla_{(\boldsymbol{\gamma}^{\top},\boldsymbol{\mu}^{\top},\boldsymbol{v}^{\top})^{\top}}f_v^*/f_v^*$ as defined in ((ref)) and

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

In ((ref)), $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)), define

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

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} ]$.

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) and (ref) in view of $LR_{n}=\lim_{\zeta\to 0}\max \{LR_{n,\zeta}^1 ,LR_{n,\zeta}^2 \}$.

This result generalizes Theorem 2 of chenchen03sinica, who derive the asymptotic distribution of the LRTS in the univariate case.

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:

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

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))--((ref)). 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

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)) but using the density ((ref)) in place of ((ref)) with the common value of $v^*$ across components; ${s}_{\alpha}^m(\boldsymbol{\lambda}_m)$ is defined as

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)) but using $\widetilde{\boldsymbol{s}}(\boldsymbol{x},\boldsymbol{z})$ defined in ((ref)) in place of ((ref)). 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

where $\Lambda_{\boldsymbol{\lambda}}^j$ is given by ((ref)).

Define $\widehat{t}_{\alpha,m}(\boldsymbol{\lambda}_m)$ by

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)), 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})$.

Simulation

Choice of penalty function

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 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)) and set $a_n = n^{-1/2}$ as recommended by 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$.

Simulation 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 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) 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). 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) reports the powers of the EM test when $a_n=1$ under three alternative models given in Table (ref). 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) 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). The EM test gives accurate type I errors across two models, sample sizes, and the values of $a_n$. Table (ref) 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). Overall, the EM test shows good power under finite sample size.

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. }

The flea beetles

The flea beetles data available in R package 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 lubischew62bio.} Figure (ref) 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), 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) 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.

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 pan02bio,he06csda. As in 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), 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 he06csda, the BIC chooses the five-component model. Table (ref) 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).

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.