On LR simultaneous test of high-dimensional mean vector and covariance matrix under non-normality

Accepted Manuscript On LR simultaneous test of high-dimensional mean vector and covariance matrix under non-normality Zhenzhen Niu, Jiang Hu, Zhidong Bai, Wei Gao

PII: DOI: Reference:

S0167-7152(18)30324-9 https://doi.org/10.1016/j.spl.2018.10.008 STAPRO 8347

To appear in:

Statistics and Probability Letters

Received date : 3 January 2018 Revised date : 2 October 2018 Accepted date : 7 October 2018 Please cite this article as: Niu Z., et al., On LR simultaneous test of high-dimensional mean vector and covariance matrix under non-normality. Statistics and Probability Letters (2018), https://doi.org/10.1016/j.spl.2018.10.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

On LR simultaneous test of high-dimensional mean vector and covariance matrix under non-normality Zhenzhen Niu, Jiang Hu∗ , Zhidong Bai, Wei Gao School of Mathematics & Statistics, Northeast Normal University, Changchun, P.R.C. Key Laboratory for Applied Statistics of Ministry of Education

Abstract In this paper, we primarily focus on simultaneous testing mean vector and covariance matrix with high-dimensional non-Gaussian data, based on the classical likelihood ratio test. Applying the central limit theorem for linear spectral statistics of sample covariance matrices, we establish new modification for the likelihood ratio test, and find that this modified test converges in distribution to normal distribution, when the dimension p tends to infinity, proportionate to the sample size n under the null hypothesis. Furthermore, we conduct a simulation study to examine the performance of the test and compare it with other tests proposed in past studies. As the simulation results show, our empirical powers are clearly superior to those of other tests in a series of settings. Keywords: high-dimension, simultaneous test, mean vector, covariance matrix, non-Gaussian distribution, RMT 2010 MSC: 62H15, 62H10

1. Introduction In the classical multivariate statistical analysis, statisticians usually consider estimation and hypothesis tests under the assumption of a large sample size n, but the dimension p of the random vector is fixed, such as Anderson (2003); 5

Muirhead (1982); Eaton (1983). However, with rapid developments in and ex∗ Corresponding

author.

Preprint submitted to Statistics & Probability Letters

October 11, 2018

tensive application of computers, the storage and analysis of high-dimensional data is possible. A common feature of all these high-dimensional data is that their dimensions can be proportionally large when compared with the sample size, such as the analysis of the stock market in economics, signal processing 10

in wireless communications, and the study of sequences of genes in biology. A special example is the simultaneous test of the mean vector and the covariance matrix, which plays an important role in multivariate statistical analysis to consider whether the sample is from some known population or not. Let {x1 , . . . , xn } be the independent and identically distributed (i.i.d.) ob-

15

servations from a p-dimensional population x with mean vector µp and covariance matrix Σp . Now, let us consider the hypothesis test H0 : µp = µ0 and Σp = Σ0

v.s.

H1 : µp 6= µ0 or Σp 6= Σ0 ,

(1.1)

where µ0 and Σ0 are a given vector and a positively defined matrix, respec−1

tively. Note that by the transformation Σ0 2 (xi − µ0 ), we can assume µ0 = 0p and Σp = Ip respectively without loss of generality. Here, and in the sequel, 0p 20

denotes the p-dimensional zero vector and Ip is the p×p identity matrix. Therefore, the above hypothesis test (1.1) is equivalent to the following hypothesis test H0 : µp = 0p and Σp = Ip

v.s.

H1 : µp 6= 0p or Σp 6= Ip .

(1.2)

If the population is normally distributed, one of the most important tests in statistics is the likelihood ratio test (LRT), which can be found in any textbook 25

on multivariate statistics analysis, for example, see Section 10.9 (page 444) in Anderson (2003). That is ¯0x ¯, N H = (trS − log |S| − p) + x

(1.3)

where n

¯= x

1X xi , n i=1

n

S=

1X ¯ )(xi − x ¯ )0 . (xi − x n i=1

Further assume that the dimension p is fixed and the sample size n goes to infinity, then, under the null hypothesis, it was proved that the classical LRT, 2

30

that is, n times N H in (1.3) converges to χ2 distribution with f = p(p+1)/2+p degrees of freedom. Moreover, Sugiura and Nagao (1968) and Das Gupta (1969) showed that the classical LRT is unbiased. These excellent properties encourage us to investigate what is going on if p is large and the population is not normally distributed. Therefore, in this paper, we focus on the high-dimensional LR

35

simultaneous test for mean vector and covariance matrix under non-normality. Before presenting our main results, we first review some existing results in the literature about testing the high-dimensional mean vector and covariance matrix. For the problem of testing µp = 0p , Pan and Zhou (2011) first obtained ¯ ¯ 0 S−1 x the central limit theorem (CLT) of the classical Hotelling’s T 2 statistic x

40

under the condition lim p/n ∈ (0, 1) and without the normality assumption. Later, Jiang and Yang (2013) and Jiang and Qi (2015) showed that under the normality assumption, the CLT still holds when p < n − 1 and p → ∞ as n goes to infinity. In addition, there are three modifications of Hotelling’s T 2

statistic by Srivastava and Du (2008), Chen and Qin (2010) and Chen et al. 45

(2011) respectively for the singular S case, which can be viewed as naive tests with Euclidean norm and can be found in Hu and Bai (2016). For one sample scatter test problem, that is, testing Σp = Ip , Bai et al. (2009), Jiang and Yang (2013) and Jiang and Qi (2015) provided the high dimensional corrections for the LRT. Ledoit and Wolf (2002) proposed a test statistic tr(S − Ip )2 and

50

showed its CLT under lim p/n ∈ (0, ∞). In order to study the simultaneous testing of high-dimensional mean vector and covariance matrix, Jiang and Yang (2013) and Jiang and Qi (2015) also considered this problem under the normality assumption. Additionally, Theorem 2 in Chen and Jiang (2018) provided the

55

CLT under an alternative hypothesis. The most recently work is by Liu et P ¯0x ¯ + tr( n1 ni=1 xi x0i − Ip )2 and al. (2017), who proposed the statistic LL = x

provided its asymptotic distribution under the null hypothesis. For more details on other test problems, we refer to the review paper Hu and Bai (2016). The rest of the paper is organized as follows. In section 2, we present the main results of the corrected LRT statistic which is defined as (2.2). In section 60

3, we show a brief simulation of the results, including the empirical sizes and 3

empirical powers. The proof of the main result as well as some useful results from the random matrix theory (RMT) are explained in Section 4. Section 5 provides some conclusions and discussions. 2. Main results 65

Before presenting our main results, we first present a few notations and some D

p

basic results found in RMT. In the following, the notation −→ and − → mean that “converge in distribution” and “converge in probability”, respectively. For any

70

n × n matrix A with only real eigenvalues, we denote F A as the empirical specPn A tral distribution (ESD) of A, that is F A (x) = n1 i=1 I(λA i ≤ x), where λi

denotes the i-th smallest eigenvalue of A and I(·) is the indicator function. For any function of a bounded variation G on the real line, its Stieltjes transform R 1 dG(λ), where z ∈ C+ = {z : =(z) > 0}. Conis defined by mG (z) = λ−z

sider a set of samples {x1 , . . . , xn } independently drawn from a p-dimensional population x with the mean vector µp , covariance matrix Σp and satisfying the

75

1/2

linear transformation model, that is, x = Σp w + µp , w = (w1 , . . . , wp )0 and {wi , i ≤ p} are i.i.d. real random variables. Let Sn =

n 1 X ¯ )(xi − x ¯ )0 , (xi − x N i=1

(2.1)

where N = n − 1, then we have the following lemma. Lemma 2.1 (Theorem 1.1 in Silverstein (1995)). Under the condition min{p, n} → ∞, the ESD of Σp , denoted by Hp , converges weakly to a distribution H and yn = p/n → y ∈ (0, ∞), the ESD of Sn almost surely converges weakly to a

distribution function F y,H with the Stieltjes transform my,H (z) satisfying the equation my,H =

Z

1 dH(λ). λ(1 − y − yzmy,H ) − z

Notice that compared with S, Sn is an unbiased estimator of Σp . Therefore, according to this property, we present a modification of LRT (1.3), which is ¯0x ¯. N H = (trSn − log |Sn | − p) + x 4

(2.2)

80

Now, we are at the right place to present our main results. Theorem 2.2. Under the assumptions that 1/2

• (Model Assumption:) xi = Σp wi + µp , for i = 1, . . . , n, where Σp is a positive definite matrix, wi = (wi1 , . . . , wip )0 and {wij , i ≤ n, j ≤ p} are i.i.d. real random variables; 85

• (Moments Assumption:) Ew11 = 0, E|w11 |2 = 1, E|w11 |4 = ∆ + 3 < ∞; • (Dimension Assumption:) yN :=

p N

→ y ∈ (0, 1) as N, p → ∞;

• (Norm Assumption:) The spectral norm of Σp is bounded in p, and the ESD of Σp (i.e., Hp ) converges weakly to a distribution H as p → ∞, and g(x) = x − log x − 1 , we have 1

ν − 2 (N H − p · F yN ,Hp (g) − 90

1 D trΣp − µ0p µp − µ) −→ N (0, 1), n

where F yN ,Hp (g) =

Z

g(x)dF yN ,Hp (x),

(2.3)

F yN ,Hp is the limit empirical distribution of Sn (i.e., F y,H ) with y, H replaced by yN , Hp respectively, R I y m3 (z)t2 (1 + tm(z))−3 dH(t) 1 R dz µ=− g(x) 2πi C [1 − y m2 (z)t2 (1 + tm(z))−2 dH(t)]2 R I y m3 (z)t2 (1 + tm(z))−3 dH(t) ∆ R − · g(x) dz, 2πi C 1 − y m2 (z)t2 (1 + tm(z))−2 dH(t)

and

I I 1 g(z1 )g(z2 ) dm(z1 )dm(z2 ) 2π 2 C1 C2 (m(z1 ) − m(z2 ))2 I I Z y∆ t − 2 g(z1 )g(z2 ) × [ 4π C1 C2 (m(z1 )t + 1)2 t × dH(t)]dm(z1 )dm(z2 ), (m(z2 )t + 1)2

(2.4)

ν =−

m =

y−1 z

(2.5)

+ ymy,H . The contours C, C1 and C2 are contained in the analytic

region g(x), and each encloses the support of the LSD F y,H , and C1 and C2 are assumed to be non-overlapping. 5

The proof of the Theorem 2.2 will be presented in section 4. 95

Remark 2.3. Notice that if y > 1, there exist large enough p and N , such that p/N > 1. In this case, the sample covariance matrix Sn almost surely has p − N zero eigenvalues. Then, the modified LRT statistic (2.2) is not well defined because of log |Sn |. It seems possible to get rid of the zero eigenvalues from log |Sn | to solve the problem for y > 1. However, we possibly need to derive

100

a new limiting theorem and thus it would be too long to include it in the present paper. We shall keep it in our further research. Corollary 2.4. Assuming the Model, Moments and Dimension Assumptions in Theorem 2.2 hold true, then under the null hypothesis (1.2), we have −1/2

νN

D

(N H − p · lN − yN − µN ) −→ N (0, 1),

where 105

N) lN = 1− yNyN−1 log(1−yN ), µN = − log(1−y + yN2∆ , νN = −2 log(1−yN )−2yN . 2

Remark 2.5. This corollary is a special case of Theorem 2.2 with µp = 0p and Σp = Ip . Notice that µ and ν in Theorem 2.2 are the limiting results which are involved with the limit of the ratio of dimension to sample size. However, for any real sample, the limit y does not exist actually. Thus, in real applications 110

one has to use yN instead of y to get estimates of µ and ν, that is, µN and νN . Therefore, we shall subscript N to indicate the distinction. Let α be the significance level. Then, the proposed test (1.2) rejects H0 −1/2

if νN

(N H − p · lN − yN − µN ) > zα , where zα is the upper α quantile of

the standard normal distribution. Specifically, the rejection region of the test 115

problem (1.2) is given by {(x1 , . . . , xn ) : N H >

√

νN · zα + p · lN + yN + µN }.

Remark 2.6. Suppose that the condition of Theorem 2.2 holds true, Then, under the alternative hypothesis H1 of the test problem (1.2), we can get the asymptotic distribution of N H under the population covariance matrix Σp and 6

the mean value µp are given, that is, −1

νA 2 (N H − p · lA − 120

1 D trΣp − µ0p µp − µA ) −→ N (0, 1), n

where lA , µA and νA are the values of F yN ,Hp (g), µ and ν under the alternative hypothesis, and can be obtained from (2.3),(2.4) and (2.5), respectively. Further, the theoretical power function is asymptotically close to √ νN · zα + p · lN + yN + µN − p · lA − n1 trΣp − µ0p µp − µA ), βL = 1 − Φ( √ νA where Φ(·) is the distribution function of the standard normal N (0, 1). In the hypothesis test (1.2), we only require that the population has finite

125

fourth moments. In practical applications, the fourth moments are unknown, and therefore, the estimation of the fourth moments is necessary. By applying the method of moments and RMT, Zhang et al. (2017) obtained the estimation of the fourth moments of the real random variables. Since the sample size n is larger than the dimension p in this paper, similar to the Theorem 2.7 in Zhang

130

ˆ et al. (2017), we obtain the following estimator ∆. Theorem 2.7. Under the same assumption of Theorem 2.1, we have that the estimator of ∆, ˆ = (1 − yN )2 ∆

Pn

j=1 [(xj

¯ )0 S−1 ¯) − −x nj (xj − x np

p 2 1−yN ]

−

2 1 − yN

(2.6)

is weakly consistent and asymptotically unbiased. Here Snj is the sample covariance matrix of the observation by removing the vector xj . 135

The proof of this theorem is analogous to Theorem 2.7 in Zhang et al. (2017), thus we omit the details here.

3. Simulations In this section, we simulate the above hypothesis test (1.2) H0 : µp = 0p and Σp = Ip . At the same time, in order to show the performance, we 140

compare the simulation results with the statistic LL proposed by Liu et al. 7

Table 1: Empirical sizes and powers of the tests N H and LL under Model I, based on 10,000 replications with real Gamma data Power under H1 =0 = NH LL NH 0.9631 0.4553 0.9613 0.5732 0.1059 0.5726 0.1463 0.0517 0.1476 0.0823 0.0577 0.0783 0.0637 0.0551 0.063 1 1 1 1 0.9349 1 0.982 0.452 0.9826

Size under H0 n=1000, p=50 n=1000, p=100 n=1000, p=300 n=1000, p=600 n=1000, p=900 n=10000, p=100 n=10000, p=300 n=10000, p=500

NH 0.0665 0.0604 0.0539 0.0525 0.049 0.0587 0.0548 0.0521

LL 0.0611 0.0537 0.0508 0.0511 0.0491 0.0577 0.0537 0.0489

0.1 LL 0.4533 0.1039 0.0511 0.0521 0.0569 1 0.9381 0.4567

(2017) under various settings. However, limited to the space of the paper, we can only present two of them herewith. What is more, Jiang and Yang (2013) has proved the poor performance of the classical LRT under normal distribution and in this part, we will not mention it any further. In the following discussion, 145

[x] represents the maximum integer that is not greater than the real number x. The random samples {xi }ni=1 are generated from the following model, that 1

is xi = µp + Σp2 wi , i = 1, . . . , n, where wi = (wi1 , . . . , wip )0 and assume that {wij , j = 1, . . . , p} are i.i.d. random variables from Gamma distribution

Gamma(4, 2)−2. Note that in the following simulations, we always assume that 150

the fourth moments of the {wij }pj=1 are known. As per the simulation results shown by Zhang et al. (2017) in Tables 5 and 6, we believe that the performance of the test is essentially invariant when the true values of the fourth moments are replaced by the estimates. In this paper, the result of the simulation between the real and the estimated values of the four moments is omitted.

155

For the mean vector µp and the covariance matrix Σp , we consider the following different models: • Model I: µp = (, 00p−1 )0 with  ∈ {0, 0.1}, and Σp = diag(0.5, 10p−1 ), where 1p denotes the p-dimensional vector with all elements being 1. • Model II: µp = 0p , and Σp = diag(1 − ρ, . . . , 1 − ρ, 1, . . . , 1), where the

160

number of (1 − ρ) is equal to [p/4] and ρ is a constant. For the two models mentioned above, we run a simulation with 10,000 rep-

8

etitions to get the empirical sizes and powers of N H and LL at a significant level of α = 0.05, where we take n ∈ {100, 1000, 10000} and choose different dimensions p. The results are showed in Table 1 and Figure 1. Table 1 presents 165

the empirical sizes and empirical powers of N H and LL under Model I. Figure 1 shows the divergence of powers of the two test statistics when the parameter ρ increases from 0 to 1 under Model II. These results show that the empirical sizes of the two tests are all close to the significance level 0.05 when the dimension p increases. From Table 1, we find that N H has far higher powers than

170

LL, indicating that N H is more sensitive than LL under Model I. Actually, we also did many other models, and the results find that when p/N is small or the eigenvalues of Σp do not have the same size, our statistic would be better. In addition, from Figure 1 we know that if p/N is close to 1, then the power of N H becomes worse than LL. That is because when p/N → 1, the LRT of

175

covariance matrix performs poorly, which can be found in Bai and Saranadasa (1996) for more discussion.

Figure 1: The empirical powers are estimated based on 10,000 replications with real Gamma variables under Model II and these results are based on the significance level of α = 0.05.

9

4. Proof of Theorem 2.2 In this section we will present a brief proof of the Theorem 2.2. We first provide a lemma, which will be useful for our proof. 180

Lemma 4.1 (Theorem 2.2 in Zheng et al. (2015)). Assume that the population x fulfills Model Assumption, Moments Assumption, Dimension Assumption, and Norm Assumption in Theorem 2.2. Let f1 , . . . , fk be functions ana-

185

lytic on an open domain of the complex plan containing the support of the LSD Pp F y,H and define Yp (fl ) = p{µSn (fl )−F yN ,Hp (fl )} = i=1 fl (λi )−pF yN ,Hp (fl ),

where {λi }pi=1 are the eigenvalues of the unbiased sample covariance matrix R Sn in (2.1) and F yN ,Hp (fl ) = fl (x)dF yN ,Hp (x). Then, the random vector (Yp (f1 ), . . . , Yp (fk )) converges to a k-dimensional Gaussian random vector

(Yf1 , . . . , Yfk ) with mean function µ and variance-covariance function ν, where µ and ν are defined in (2.4) and (2.5) respectively. ¯. Proof of Theorem 2.2. Recall the LRT statistic N H = (trSn −log |Sn |−p)+¯ x0 x By the definition of ESD of Sn , denoted by Fn in the following, Z ¯0x ¯ N H = p (x − log x − 1)dFn (x) + x Z ¯0x ¯ = p g(x)d(Fn (x) − F yN ,Hp (x)) + p · F yN ,Hp (g) + x

190

g ¯0x ¯, =N H +x

g where N H =p

R

g(x)d(Fn (x) − F yN ,Hp (x)) + p · F yN ,Hp (g). It is easy to see

that the proof of Theorem 2.2 can be divided into two parts, one is to prove R p g(x)d(Fn (x)−F yN ,Hp (x)) converges in distribution to a normal distribution,

¯0x ¯ converges in probability to a constant. The and the other is to prove that x

195

g result in Lemma 4.1 indicates that N H − p · F yN ,Hp (g) converges weakly to a D g Gaussian distribution, that is N H − p · F yN ,Hp (g) −→ N (µ, ν), where the mean

µ and variance function ν are defined in (2.4) and (2.5), respectively. Thus, the p

¯0x ¯ − n1 trΣp − µ0p µp − → 0, in which we only need remainder serves to prove that x

to prove ¯) = E(¯ x0 x

1 trΣp + µ0p µp , n 10

¯ ) → 0. Var(¯ x0 x

(4.1)

In fact, (4.1) has been proved in the process of proving Theorems 2.1 and 2.2 200

in the Appendix A of Liu et al. (2017), when µ0p Σp µp tends to a constant and n−1 trΣ2p converges, and thus we omit the proof here. By summarizing the above results, and using Slutsky lemma, we complete the proof of Theorem 2.2.

5. Conclusions and discussions In this paper, we propose corrected LRT for the simultaneous testing of 205

mean vector and covariance matrix for high-dimensional non-Gaussian data. Motivated by past studies by Jiang and Yang (2013), who provide the CLTs for several classical LR tests for mean vectors and covariance matrices of highdimensional Gaussian distributions, we provide the asymptotic distribution of the LRT under the null hypothesis and provide the theoretical power under

210

the alternative hypothesis, by applying the method of RMT. Simulation results show that our LRT performs well. It is noteworthy that this paper only studies the test of a single population, and the simultaneous testing of the mean vectors and covariance matrices of multiple populations will be our focus in the future.

Acknowledgment 215

This work was partially supported by CNSF 11771073, 11571067, and 11471068. The authors would like to thank the two anonymous referees, whose comments improve the quality of the paper significantly.

References Anderson, T. (2003). An introduction to multivariate statistical analysis. Third 220

Edition. New York: Wiley. Bai, Z. D. and Saranadasa, H. (1996). Effect of high dimension: by an example of a two sample problem. Statistica Sinica, 311-329.

11

Bai, Z. D., Jiang, D. D., Yao, J. F. and Zheng, S. R. (2009). Corrections to LRT on large-dimensional covariance matrix by RMT. The Annals of Statistics, 225

37(6B), 3822-3840. Chen, H. J. and Jiang, T. F. (2018). A study of two high-dimensional likelihood ratio tests under alternative hypotheses. Random Matrices: Theory and Applications, 7(01), 1750016. Chen, S. X. and Qin, Y. L. (2010). A two-sample test for high-dimensional

230

data with applications to gene-set testing. The Annals of Statistics, 38(2), 808-835. Chen, L. S. Paul, D., Prentice, R. L. and Wang, P. (2011). A regularized Hotelling’s T 2 test for pathway analysis in proteomic studies. Journal of the American Statistical Association, 106(496), 1345-1360.

235

Das Gupta, S. (1969). Properties of power functions of some tests concerning dispersion matrices of multivariate normal distributions. The Annals of Mathematical Statistics, 40(2), 697-701. Eaton, M. L. (1983). Multivariate statistics: A vector space approach[M]. Wiley, New York.

240

Hu, J. and Bai, Z. D. (2016). A review of 20 years of naive tests of significance for high-dimensional mean vectors and covariance matrices. Science China Mathematics. 59(12), 2281-2300. Jiang, T. F. and Qi, Y. C. (2015). Likelihood ratio tests for high-dimensional normal distributions. Scandinavian Journal of Statistics, 42(4), 988-1009.

245

Jiang, T. F. and Yang, F. (2013). Central limit theorems for classical likelihood ratio tests for high-dimensional normal distributions. The Annals of Statistics, 41(4), 2029-2074. Ledoit, O. and Wolf, M.(2002). Some hypothesis tests for the covariance matrix when the dimension is large compared to the sample size. The Annals of

250

Statistics, 30, 1081-1102. 12

Liu, Z. Y., Liu, B. S., Zheng, S. R. and Shi, N. Z. (2017). Simultaneous testing of mean vector and covariance matrix for high-dimensional data. Journal of Statistical Planning and Inference, 188, 82-93. Muirhead, R. J. (1982). Aspects of Multivariate Statistical Theory. Wiley, New 255

York. Pan, G. M. and Zhou, W. (2011). Central limit theorem for Hotelling’s T 2 statistic under large dimension. The Annals of Applied Probability, 21(5), 1860-1910. Silverstein, J. W. (1995). Strong convergence of the empirical distribution of

260

eigenvalues of large dimensional random matrices, Journal of Multivariate Analysis, 55(2), 331-339. Srivastava, M. S. and Du, M. (2008). A test for the mean vector with fewer observations than the dimension. Journal of Multivariate Analysis, 99(3), 386-402.

265

Sugiura, N. and Nagao, H. (1968). Unbiasedness of some test criteria for the equality of one or two covariance matrices. The Annals of Mathematical Statistics, 39(5), 1686-1692. Zhang, Q. Y., Hu, J. and Bai, Z. D. (2017). Optimal modification of the LRT for the equality of two high-dimensional covariance matrices. arXiv: 1706.

270

06774v3. Zheng, S. R., Bai, Z. D. and Yao, J. F. (2015). Substitution principle for CLT of linear spectral statistics of high-dimensional sample covariance matrices with applications to hypothesis testing. The Annals of Statistics, 43(2), 546-591.

13