Testing the equality of a large number of means of functional data

Journal of Multivariate Analysis 185 (2021) 104778 Contents lists available at ScienceDirect Journal of Multivariate Analysis journal homepage: www.elsevier.com/locate/jmva Testing the equality of a large number of means of functional data M. Dolores Jiménez-Gameroa, Alba M. Franco-Pereirab,∗ aDpto. Estadística e Investigación Operativa, Universidad de Sevilla, Spain bDpto. Estadística e Investigación Operativa, Universidad Complutense de Madrid, Spain article info Article history: Received 23 July 2020 Received in revised form 18 May 2021 Accepted 18 May 2021 Available online 25 May 2021 AMS 2010 subject classifications: 62G10 Keywords: Asymptotic normality Consistency Finite dimensional data Functional data Testing the equality of means abstract Given kindependent samples of functional data, this paper deals with the problem of testing for the equality of their mean functions. In contrast to the classical setting, where kis kept fixed and the sample size from each population increases without bound, here kis assumed to be large and the size of each sample is either bounded or small in comparison to k. A new test is proposed. The asymptotic distribution of the test statistic is stated under the null hypothesis of equality of the kmean functions as well as under alternatives, which allows us to study the consistency of the test. Specifically, it is shown that the test statistic is asymptotically free distributed under the null hypothesis. The finite sample performance of the test based on the asymptotic null distribution is studied via simulation. Although we start by assuming that the data are functions, the proposed test can also be applied to finite dimensional data. The practical behavior of the test for one dimensional data is numerically studied and compared with other tests. ©2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Let X1,...,Xkbe krandom vectors taking values in Rd, for some d∈N, with finite means and variances. Let Xi,1,...,Xi,nibe i.i.d. (independent and identically distributed) from Xi, 1 ≤i≤k. The ksamples are assumed to be independent. A classical problem in Statistics is that of testing for the equality of the means of the kpopulations, that is, testing for H0:µ1= ··· = µk, with µj=E(Xj), 1 ≤j≤k. When d=1, if the populations are homoscedastic, following normal distributions and kis a fixed quantity, then the ANOVA F-test provides an exact solution to the above problem. If the normality assumption is dropped but the sample sizes are large, then the F-test is asymptotically valid. For heteroscedastic populations, Cochran’s test gives an asymptotically free distributed test. A similar scenario could be described for d-variate data. In the context of functional data, it is also of interest testing for the equality of mean functions. Some papers have dealt with this problem: Benko et al. [4] (for k=2), Horváth and Kokoszka [19] (for k=2) and Horváth and Rice [20] (for k≥2) have proposed tests based on projections onto the space determined by the leading eigenfunctions of an estimator of the covariance operator of the joint population; Ghiglietti et al. [11] (for k=2) have proposed a test based on a generalized Mahalanobis distance for infinite dimensional spaces whose computation involves calculating the eigenvalues and eigenfunctions of an estimator of the covariance operator of the joint population; Yuan et al. [28] and Zhang et al. [33] ∗Corresponding author. E-mail address: [email protected] (A.M. Franco-Pereira). https://doi.org/10.1016/j.jmva.2021.104778 0047-259X/©2021 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/). M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 (for k=2), and Cuevas et al. [9], Zhang et al. [31,32] and the approaches reviewed in Zhang [30] (for k≥2) consider tests based on calculating the norms of the difference between sample means; Cuesta-Albertos and Febrero-Bande [8] (for k≥2) have proposed a test based on random projections; among many others. A common way of approximating the null distribution of a test statistics is by means of its asymptotic null distribution, whenever it does not depend on unknown quantities. The test statistics in [4,19,20] are (under the null) asymptotically free distributed, although the associated tests are not universally consistent. Other test statistics have asymptotic null distributions depending on unknowns and, hence, their null distribution must be approximated by using alternative estimators. The consistency of the bootstrap approximation has been stated in Paparoditis and Sapatinas [24]. So far we have assumed that kis arbitrary but fixed and used the term asymptotically to mean ni→ ∞, 1 ≤i≤k. The case in which kis allowed to increase has also received attention in the statistical literature. In particular, when d=1, if the sample sizes n1= ··· = nkare fixed and assuming that the data are i.i.d., Boos and Brownie [5], and Akritas and Arnold [1] have proven that, conveniently normalized, the F-test has asymptotically a standard normal distribution. This result is also true for unbalanced designs, but in such a case the asymptotic variance depends on the kurtosis of the population (the fourth moment is assumed to exist). Harrar and Gupta [17] devised a transformation that speeds up the convergence to the limiting distribution. Assuming independent data from a d-variate population (d≥2) with the same distribution across the groups except for their means, which are allowed to vary between groups, Gupta et al. [15] have proven that, conveniently normalized, commonly used statistics in MANOVA are also asymptotically normal distributed. Besides, they proved that the variance of the asymptotic distribution depends on fourth order moments of the population, which are assumed to exist. Although it is not required the sample sizes to be equal, they are assumed to be close (in certain sense). Notice that those papers assume, among other things, that the populations have equal variances. For univariate normal populations, that may have different variances, when the samples have ‘‘similar’’ sample sizes (in the sense that maxini<aminini, for some positive constant a), Park and Park [25] proposed two tests whose associated statistics, conveniently normalized, are asymptotically normal. The assumption that the populations have all of them a normal distribution is crucial in order to derive the asymptotic distribution. In fact, we have carried out some simulations, with data coming from a negative exponential distribution (and other non-normal populations) for several values of the common mean and several values of k, and observed that the empirical type I errors were far apart from the nominal value (results are reported in Section 7). Akritas and Papadatos [2] and Harrar and Bathke [16] proposed tests for the equality of a large number of univariate and d-variate means, not requiring homoscedasticity, respectively. Assumption 1 in [16], needed to derive the asymptotic null distribution of the tests statistics in that paper, and some of these test statistics cannot be carried to the infinite dimensional case. This paper proposes and studies a test of H0with klarge, requiring neither normality nor homoscedasticity nor equality of distributions. We start by assuming that the data are functions. Our procedure is fully functional, not based on projections. The test statistic is asymptotically free distributed, not relying on resampling or Monte-Carlo methods to obtain critical values. Then, we show the way the new test can be applied to data with finite (but arbitrary) dimension. The sample sizes can be bounded or they can increase with kat a certain rate that will be specified later in the text. The design can be balanced or unbalanced but, for unbalanced designs, the ratio between the maximum and the minimum sample size cannot be greater than a certain power of k. One of the advantages of our proposal is that, in comparison with existing tests for finite dimensional data, it requires weaker assumptions. The paper unfolds as follows. Based on a characterization of H0, Section 2proposes a test statistic. In Section 3it is shown that, adequately normalized, it asymptotically has a standard normal distribution. The practical use of such statistic requires an estimator of its null variance. This issue is dealt with in Section 4. The power of the test rejecting the null for large values of the test statistic is studied in Section 5, in which alternatives that this test is able to detect are identified. The results in Sections 2–5 are valid for k→ ∞. Section 6summarizes the results of a simulation study designed to study the goodness of the asymptotic null distribution as an approximation of the null distribution of the test statistic, and to study the power of the proposed test, for a finite k. This section also contains two real data set applications. In the previous paragraph we cited some tests for the comparison of the means of k(≥2) independent populations, that were designed for kfixed and increasing sample sizes. One may wonder if such tests are still valid in the setting considered in this paper. Section 7deals with this issue for the tests proposed in [8,9,20,31,32] and the approaches reviewed in [30]. This section also investigates the validity of the new test in the case of a fixed kand large sample sizes. Although in those sections we have assumed functional data, Section 8shows how the test can be applied to finite dimensional data. This section also displays the results of some simulation experiments, where the proposed test is compared with those in [25] for one-dimensional data. A real data set application is included as well. Section 9concludes the paper. All proofs are deferred to Section 10. All computations have been carried out using programs written in the Rlanguage [26]. The Rcode for the calculation of the proposed test statistic is available in the Supplementary Material. The following notation will be used along the paper: all limits in this paper are taken when k→ ∞;L →denotes convergence in distribution; P →denotes convergence in probability; a.s. →denotes the almost sure convergence; Mdenotes a positive constant whose exact value is unimportant and may vary across the text; and an unspecified integral denotes integration over the compact interval [0,1]. 2 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 2. The test statistic Assume that X1,...,Xkare random functions defined on a common probability space (Ω,A,P) with values in the separable Hilbert space L2=L2([0,1],R), the space of square-integrable real-valued functions defined on the compact interval [0,1], with the usual inner product ⟨f,g⟩ = ∫f(t)g(t)dt and norm ∥f∥ = ⟨f,f⟩1/2,f,g∈L2. Assume that Xiis integrable, that is, that E∥Xi∥<∞, which implies that there is a unique function µi∈L2, the mean function of Xi, satisfying E(⟨Xi,x⟩)= ⟨µi,x⟩,∀x∈L2, 1 ≤i≤k. It follows that µi(t)=E{Xi(t)}for almost all t∈ [0,1]. Let ci(t,s)=E[{Xi(t)−µi(t)}{Xi(s)−µi(s)}],s,t∈ [0,1], stand for the covariance function of Xi, and let Ci:L2→L2denote the covariance operator of Xi, defined as Cif(t)=∫ci(t,s)f(s)ds,∀f∈L2, which is assumed to satisfy 0< θi=∫∫ c2 i(t,s)dtds <∞,1≤i≤k.(1) Along this paper it is assumed that the available data consist of independent random samples from each population, that is, we assume that Xi= {Xi1,...,Xini}are niindependent random functions from Xi, 1 ≤i≤k, and that X1,...,Xk are independent. On the basis of the available data, we are interested in testing H0, where the equalities under H0are understood in the L2sense, that is, ∥µi−µj∥ = 0, ∀i= j. Notice that H0is equivalent to Dk=0, where Dk=1 k k ∑ i=1∥µi− ¯µ.∥2=1 k k ∑ i=1⟨µi, µi⟩−⟨¯µ.,¯µ.⟩,¯µ.=1 k k ∑ i=1 µi. Since the sample mean is unbiased for the population mean, the above equivalence leads us to consider the following test statistic for testing H0, ˜ Tk=1 k k ∑ i=1∥¯ Xi.−¯ X..∥2,with ¯ Xi.=1 ni ni ∑ j=1 Xij,1≤i≤k,¯ X.. =1 k k ∑ i=1¯ Xi.. Routine calculations show that ˜ Tk=(1 k−1 k2)k ∑ i=1˜ h1(Xi)−1 k2 k ∑ i,l=1 i=l h2(Xi,Xl), with ˜ h1(Xi)=h2(Xi,Xi), h2(Xi,Xl)= ⟨¯ Xi.,¯ Xl.⟩. We have that E{h2(Xi,Xl)}=⟨µi, µl⟩,∀i= l, and, assuming ∫ci(t,t)dt <∞, E{˜ h1(Xi)}=⟨µi, µi⟩+ 1 ni∫ci(t,t)dt. Therefore ˜ Tkis a biased estimator of Dk. The bias can be removed by replacing ˜ h1(Xi) with h1(Xi), defined as h1(Xi)=1 ni(ni−1) ni ∑ u,v=1 u=v ⟨Xiu,Xiv⟩. Let Tk=(1 k−1 k2)k ∑ i=1 h1(Xi)−1 k2 k ∑ i,l=1 i=l h2(Xi,Xl). Then, E(Tk)=Dk. Since E(Tk)=0 under H0and E(Tk)>0 under alternatives, the null hypothesis is rejected for large values of Tk. Now, in order to determine what large values are in this context, we have to calculate its distribution under the null hypothesis, or (at least) an approximation of it. The null distribution of Tkis clearly unknown. Next, we try to approximate it by means of its asymptotic null distribution. 3. The asymptotic distribution of the test statistic By applying the Decomposition Lemma of Efron and Stein [10], after some rearrangements, one gets that Tk=Dk+Tk,Lin +Rk,(2) with Tk,Lin =(1/k)∑k i=1Lik,Lik =h1(Xi)−⟨µi, µi⟩−2⟨¯ Xi.−µi,¯µ.⟩, 1 ≤i≤k, and Rk= − 1 k2 k ∑ i=1{h1(Xi)−⟨µi, µi⟩}+2 k2 k ∑ i=1⟨¯ Xi.−µi, µi⟩− 1 k2 k ∑ i,l=1 i=l ⟨¯ Xi.−µi,¯ Xl.−µl⟩. Clearly, E(Tk,Lin)=E(Rk)=0. The variance of Tk,Lin is given in the next lemma. 3 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 Lemma 1. Suppose that X1,...,Xkare integrable satisfying (1), that Xi= {Xi1,...,Xini}are ni≥2independent random functions from Xi,1≤i≤k, and that X1,...,Xkare independent. Then var(Tk,Lin)=1 k2 k ∑ i=1(2 ni(ni−1)θi+4 ni γi), where γi=∫∫ ci(t,s){µi(t)− ¯µ.(t)}{µi(s)− ¯µ.(s)}dtds,1≤i≤k. Notice that var(Tk,Lin)>0 is equivalent to ∑k i=1θi>0, and the latter is implied by (1). If H0is true, then var(Tk,Lin) simplifies to var0(Tk,Lin)=2 k2 k ∑ i=1 1 ni(ni−1)θi.(3) The next result shows that the term Rkin (2) is negligible in comparison to Tk,Lin. Lemma 2. Suppose that the assumptions in Lemma 1are fulfilled. Then var(Rk)/var(Tk,Lin)→0. Since Tk,Lin is an average of independent random variables it readily follows that, conveniently scaled and under some assumptions (specifically, that the Lindeberg condition in (4) holds true), it is asymptotically normally distributed. This fact together with the result in Lemma 2 leads to the following result. Theorem 1. Suppose that the assumptions in Lemma 1are fulfilled and that 1/k2 var(Tk,Lin) k ∑ i=1 E[L2 ikI{L2 ik > ε k2var(Tk,Lin)}] → 0,∀ε > 0.(4) Then Tk−Dk √var(Tk,Lin) L −→ Z,(5) where Z has a standard normal distribution. As an immediate consequence of Theorem 1, taking into account that H0is equivalent to Dk=0, we next derive the asymptotic null distribution of Tk. Corollary 1. Suppose that the assumptions in Theorem 1are fulfilled and that H0is true. Then √kTk/σ0k L −→ Z, where Z has a standard normal distribution and σ2 0k=kvar0(Tk,Lin), with var0(Tk,Lin)as defined in (3). Recall that the test is one-sided rejecting the null hypothesis for large values of Tk. If σ2 0kwere a known quantity, in view of Corollary 1, the test that rejects H0when √kTk/σ0k≥z1−α, for some α∈(0,1), where Φ(z1−α)=1−αand Φstands for the cumulative distribution function of the standard normal distribution, would have (asymptotic) level α. In general, checking that Lindeberg condition holds is not an easy task. Because of this reason, in most cases it is easier to see that some sufficient conditions are met. Next we provide some conditions for (4) to hold. First we list such conditions, then we comment on them, and finally Proposition 1 shows that they imply (4). Assumption 1. ∃τ > 0 and k0=k0(τ)∈Nsuch that 1 k∑k i=1θi≥τ,∀k≥k0. Assumption 2. E(∥Xi−µi∥4)≤M,∀i. Assumption 3. n4 max n4 min 1 k=o(1), where nmax =max 1≤i≤kniand nmin =min 1≤i≤kni. Assumption 4. Either γ1= ··· = γk=0, or not all γisare equal to 0 and (1/k)Γk=o(1), where Γk= (1/k){∑k i=1E(⟨Xi−µi, µi− ¯µ.⟩4)/n2 i}/(1 k∑k i=1γi/ni)2 . 4 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 Assumption 1 is tantamount to say that there is a percentage of 100 ·p% of covariance operators, for some 0 <p≤1, say θ1, . . . , θ[pk], such that θi≥c, for some c>0, 1 ≤i≤ [nk]. This is clearly a nonrestrictive assumption. Assumption 2 is similar or even weaker than other comparable conditions used in the finite dimensional case: for univariate data, if the design is not balanced, the derivations in Boos and Brownie [5] require the fourth moment to exist; for data with dimension greater than 1 but finite, Assumption 1 in Gupta et al. [15] is stronger since it requires the finiteness of the moment of order 4 +δof the absolute value of each component of the random vector, for some δ > 0. Assumption 3 tells us how different the sample sizes can be in relation to k. In particular, for either balanced designs, n1= ··· = nk, or bounded sample sizes, it imposes no restriction. In the finite dimensional setting, for univariate normal populations, Park and Park [25] assume that maxini<aminini, for some positive constant a, which is more limiting than our assumption. Finally, we will see that Assumption 4 is not a severe restriction. Under H0it is automatically satisfied because γ1= ··· = γk=0. Thus, it only happens to be a constraint under alternatives. Let Yi= ⟨Xi−µi, µi− ¯µ.⟩/√ni1≤i≤k. Notice that Y1,...,Ykare kindependent random variables. Let Wkbe an equal mixture of the random variables Y1,...,Yk. Then Γkis the kurtosis of Wk. This assumption allows the kurtosis of Wkto increase with k, but at a lower rate than k, which seems not to be a hard limitation. Proposition 1. Suppose that the assumptions in Lemma 1and Assumptions 1–4are fulfilled. Then (4) holds. 4. Estimating the null variance In order to get a critical region for testing H0using the result in Corollary 1, we need a consistent estimator of σ0k. Recall that σ2 0k=k×var0(Tk,Lin)=2/k∑k i=1θi/ni(ni−1). σ2 0kis unknown because θ1, . . . , θkare unknown quantities. From expression (1),θican be unbiasedly estimated by replacing c2 i(t,s) with an unbiased estimator. Assume that ni≥4, 1≤i≤k, then c2 i(t,s) can be unbiasedly estimated by ˆ c2 i(t,s)=1 ni(ni−1)(ni−2)(ni−3) ∑ 1≤u=v=w=z≤ni h(Xiu,Xiv;t,s)h(Xiw,Xiz;t,s), where h(Xiu,Xiv;t,s)=1 2{Xiu(t)−Xiv(t)}{Xiu(s)−Xiv(s)}. Note that ˆ c2 i(t,s)=1 ni(ni−1) ∑ 1≤u=v≤ni h(Xiu,Xiv;t,s)S2 i(u,v)(t,s),(6) where S2 i(u,v)(t,s)=1 ni−3∑w=u,v{Xiw(t)−¯ Xi(u,v)(t)}{Xiw(s)−¯ Xi(u,v)(s)}, and ¯ Xi(u,v)(t)=1 ni−2∑w=u,v Xiw(t). Formula (6) is useful for computational purposes. Let ˆσ2 0k=(2/k)∑k i=1ˆ θi/ni(ni−1), where ˆ θi=∫∫ˆ c2 i(t,s)dtds. Next proposition shows that, under some mild assumptions, the quotient ˆσ2 0k/σ2 0kconverges in probability to 1. Proposition 2. Suppose that the assumptions in Lemma 1and Assumptions 1–3are fulfilled, and that ni≥4,1≤i≤k. Then ˆσ2 0k/σ2 0k P −→ 1. As an immediate consequence of Corollary 1and Proposition 2, we have the following result. Theorem 2. Suppose that the assumptions in Corollary 1and Proposition 2are fulfilled. Then √kTk/ˆσ0k L −→ Z, where Z has a standard normal distribution. Let α∈(0,1). As a consequence of Theorem 2, the test rejecting H0when √kTk/ˆσ0k≥z1−αis asymptotically correct in the sense that its type I error is asymptotically equal to the nominal value α. Remark 1. Let R(X11,...,X1n1,...,Xk1,...,Xknk)=√kTk/ˆσ0k. Notice that R(aX11 +µ, . . . , aX1n1+µ, . . . , aXk1+ µ, . . . , aXknk+µ)=R(X11,...,X1n1,...,Xk1,...,Xknk), ∀a∈R,a= 0, ∀µ∈L2, that is, the test statistic is invariant under location and scale transformations. Before ending this section, we want to underline that ˆσ0kis consistent for σ0kunder quite weak assumptions. If stronger conditions are assumed, then other estimators of σ0kmay be considered. For instance, if we could reasonably assume that X1,...,Xkhave a common variance function, c1(·,·)= ··· = ck(·,·):= c(·,·), then σ2 0k=k×var0(Tk,Lin)= (θ/k)∑k i=11/ni(ni−1), where θ=∫∫c2(t,s)dtds. In this case we have to estimate θ. Under H0, if E[{Xi(t)−µi(t)}2{Xi(s)− µi(s)}2] ≤ M, 1 ≤i≤k, 0 ≤t,s≤1, for some positive constant M, and if ∫∫c(t,t)dt <∞, then θcan be consistently estimated by means of ˆ θ=∫∫ˆ c(t,s)2dtds, where ˆ c(t,s)=(N−k)−1∑ i,j{Xij(t)−¯ Xi.(t)}{Xij(s)−¯ Xi.(s)}. A sketch of the proof of the consistency of ˆ θcan be found in Section 10. 5 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 5. Power This section studies what sort of alternatives can be detected by the proposed test. Let σ2 k=k×var(Tk,Lin). Notice that 0≤σ2 0k/σ2 k≤1.(7) Under the assumptions in Theorem 1 and Proposition 2, we have that P(√kTk ˆσ0k >z1−α)=P(√kTk−Dk σk >ˆσ0k σk z1−α−√kDk σk)≈Φ(√kDk σk−σ0k σk z1−α).(8) In the light of (7) and (8), and taking into account that Dk≥0, we consider the following three cases: (i) √kDk/σk→0, (ii) √kDk/σk→δ∈(0,∞), and (iii) √kDk/σk→ ∞. In case (iii) it is clear that P(√kTk/ˆσ0k>z1−α)→1, and thus the test is consistent against that sort of alternatives. From (8), it is apparent that to derive the power in cases (i) and (ii) we must first study the quotient σ0k/σk. Next proposition shows that, in these cases, σ0k/σk→1. Proposition 3. Suppose that the assumptions in Lemma 1and Assumptions 1and 2are fulfilled, that n2 max =o(k),(9) and √kDk/σk→δ∈ [0,∞). Then σ0k/σk→1. As an immediate consequence of Theorem 1 and Propositions 2 and 3, we have the following. Corollary 2. Suppose that the assumptions in Theorem 1and Propositions 2and 3are fulfilled. Then √kTk/ˆσ0k L −→ Z+δ, where Z has a standard normal distribution. From Corollary 2, in cases (i) and (ii) we have that P(√kTk ˆσ0k >z1−α)→Φ(δ−z1−α)={αif δ=0, > α if δ∈(0,∞). Summarizing, the proposed test asymptotically detects those alternatives such that √kDk/σk→δ∈(0,∞]. It is consistent against those alternatives satisfying √kDk/σk→ ∞. The alternatives fulfilling √kDk/σk→δ∈(0,∞) play the role of contiguous alternatives in the classical setting of fixed kand large sample sizes. The test is not able to detect those cases in which √kDk/σk→0. This is not surprising since it entails that the between means variance Dkis much more smaller than the standard deviation of Tk. As an illustration, we consider the case where 1 ≤m=mk≤kpopulations have an equal mean, say µ= 0, and the other k−mpopulations have a mean equal to 0. Let pk=m/k. For simplicity, we will also assume that the populations are homoscedastic, that is, ci(t,s)=c(t,s), 1 ≤i≤k,∀t,s∈ [0,1], and that the design is balanced, that is, ni=n, 1≤i≤k. Let γ=∫∫µ(t)µ(s)c(t,s)dtds,θ=∫∫c2(t,s)dtds, and assume that both quantities are positive. With this notation, δk:= √kDk/σk=√kpk(1 −pk)∥µ∥2/√2θ/n(n−1) +4γpk(1 −pk)/n. In this setting, the question that naturally arises is how small (or big) must mbe so that the test can detect that H0is not true? In order to answer this question, we will assume that condition (9) in Proposition 3 holds. If mis a fixed number that does not depend on k(or equivalently k−mis a fixed number that does not depend on k), then δk→0. Therefore, the test it is not able to detect alternatives with a finite number of different means. This is a special case of pk→0. Next, we see that the test can detect alternatives with pk→0 but at a lower speed than 1/k. The following cases can arise: (a) pk→p∈(0,1): then δk→ ∞, for any choice of n, and thus we are in case (iii); (b) pk→0 and nremains bounded: then the asymptotic behavior of δkcoincides with that of √kpk, and thus, (b.1) if √kpk→0, which is equivalent to pk=o(k−1/2), then we are in case (i), (b.2) if √kpk→δ∈(0,∞), then we are in case (ii), (b.3) if √kpk→ ∞, then we are in case (iii); (c) pk→0, n→ ∞,npk→ρ∈(0,+∞]: then δk→ ∞ and thus we are in case (iii). (d) pk→0, n→ ∞,npk→0: then the asymptotic behavior of δkcoincides with that of √knpk, which from (9) is √knpk=o(kpk), and thus, (d.1) if √knpk→0, then we are in case (i), (d.2) if √knpk→δ∈(0,∞), then we are in case (ii), (d.3) if √knpk→ ∞, then we are in case (iii). 6 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 The answer to the question how small must m be so that the test can detect that H0is not true is that we must take m=O(f(k)) and n=O(√k/f(k)), with fsuch that f(k)→ ∞,f(k)/k→0 and n→ ∞. For instance, we could take m=kaand n=k0.5−a, for some 0 <a<0.5, or m=log(k) and n=√k/log(k), among many possible choices. Remark 2. With the aim of improving the asymptotic power of the test in case (d.2), an anonymous reviewer suggested us to consider a max-type statistic, instead of Tk. Reasoning as in Section 2,H0is equivalent to maxi=j∥µi−µj∥2=0. Now, since ∥µi−µj∥2can be unbiasedly estimated by means of Tij =h1(Xi)+h1(Xj)−2h2(Xi,Xj), one could take Tmax =maxi=jTij as the test statistic. Let Tall =(T12,T13,...,Tk−1k)⊤∈RD, where the superindex ⊤denotes transpose and D=k(k−1)/2. If ni/n1→τi∈(0,∞)∀i, when minini→ ∞, then for each fixed k,√n1Tall converges in law to a D-variate normal random vector, Z=(Z1,...,ZD)⊤, with mean µall =(∥µ1−µ2∥2,...,∥µk−1−µk∥2)⊤and a certain covariance matrix. Therefore, Tmax converges in law to maxiZi. To derive the asymptotic distribution of Tmax we have assumed that minini→ ∞ (not required to get the asymptotic distribution of Tk). In addition, the law of maxiZiis not free of nuisance parameters, even under the null hypothesis, and thus cannot be used to approximate the null distribution of Tmax. A sound study of Tmax deserves to be considered in a separate paper. Remark 3. All our results have been stated under the tacit assumption that realizations of X11,...,Xknk, i.e., complete trajectories of functions are observable. In practice, these functions are observed at a finite grid of points and the curves X11,...,Xknkare recovered by using nonparametric techniques, such as local linear regression. The statistic is then calculated from ˆ X11,...,ˆ Xknk, which stand for the resulting curve estimators. Under suitable assumptions, all previous results remain valid when the test statistic is calculated from ˆ X11,...,ˆ Xknk. See, e.g., Jiang et al. [21], in particular the comments made after the proof of their Theorem 2. 6. Simulation results for a finite number of populations Recall that the proposed test, which rejects H0when √kTk/ˆσ0k≥z1−α, has asymptotically type I error equal to α. To study the level of the proposed test for a finite value of k, we have carried out a simulation experiment as follows. We have generated a sample from each of kpopulations (k=30,40,50,100,200,300,400,500), all of them with the same mean that, without loss of generality since the test statistic is location invariant, we took equal to 0. The size of each sample has been randomly generated from a discrete uniform random law UD{a,a+1,...,b}, with (a,b)=(5,10), (11,20). The data {Xij(t),t∈ [0,1],1≤j≤ni,1≤i≤k}have been generated in discretized versions Xij(tr), for r∈ {1,...,50}, where the values trwere chosen equispaced in the interval [0,1]. Four different types of random functions have been considered: the data were generated from a Wiener process with dispersion parameter σ=1 for all populations (recall from Remark 1 that the test statistic is invariant under scale transformations, so the same results would be obtained for other choices of σ), and three heteroscedastic versions, where the data from population iwere generated from a Wiener process with standard deviation σi, which in turn was randomly generated from a uniform distribution on (1,s), 1 ≤i≤k, s=1.5,2,3. Each case was run 10,000 times. Table 1 shows the fractions of p-values less than or equal to 0.05 and 0.10, which are the estimated type I error probabilities for nominal significance level α=0.05 and 0.10, respectively. Looking at this table we observe that, in general, the actual levels are not far from the nominal values even for small values of k, specially in the case of homoscedastic populations; as the value of σis allowed to vary in a wider interval, a larger value of kis required for the empirical levels to closely match the nominal significance values. We also see that the approximation is better for larger sample sizes. To study the power, we have considered the same type of random functions as those used for the level, but now the mean function for population iwas taken µi(t)=i kI(t>t0), 1 ≤i≤k, for several values of t0. Although the test statistic is invariant under scale and location changes, the power is greatly influenced by the choice of the scale. Recall that the power is (asymptotically) an increasing function of √kDk/σk, and thus the power increases with Dk(which measures how different the means are), with k, with the sample sizes (recall that σkdecreases as the sample sizes increase), and decreases when the scale increases. Because of this reason, in order to study the power we have considered several values of the scale parameter for the homoscedastic case. Table 2 shows the fractions of p-values less than or equal to 0.05 and 0.10, which are the estimated powers for nominal significance level α=0.05 and 0.10, respectively. They have been obtained by generating 2,000 samples in each case (for the power we observed, in some preliminary simulations, that the results using 2,000 samples and those obtained with 10,000 samples were very close, so in order to save computation time, we took 2,000 samples). Looking at this table we observe that the above considerations, which are valid asymptotically, also hold for finite sample sizes. We illustrate the applicability of our methodology through two real data examples. First, we have applied the new test to compare the mean of the occupancy rate across the 52 Spanish provinces. The data were taken from the website of the Spanish National Institute of Statistics, http://www.ine.es. For each province, we took as sample the monthly observed occupancy rate from 2011 to 2019. Therefore, we have k=52 populations, each sample curve is observed in 12 points (months) and the sample sizes are n1=n2= ··· = n52 =9. We can consider the data as i.i.d. since the accompanying note states ‘‘due to different updates in the Establishments Directory, data from different years are not directly comparable’’. Fig. 1 displays the sample curves of the occupancy rate for each province and the group sample means. Each color 7 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 Table 1 Observed proportion of rejections in 10,000 simulated data sets for each scenario, for functional data, under the null hypothesis. The size of each sample has been randomly generated from a discrete uniform random law UD{a,a+1,...,b}, with (a,b)=(5,10), (11,20). The data were generated from a Wiener process with mean 0 and with dispersion parameter σ=1 for all populations and three heteroscedastic versions. The nominal significance levels are α=0.05 and 0.10. Number of populations k Sample sizes limits ni Dispersion parameters of the Wiener processes σi=1σi∼U(1,1.5) σi∼U(1,2) σi∼U(1,3) 5% 10% 5% 10% 5% 10% 5% 10% 30 5–10 .0717 .1176 .0714 .1143 .0753 .1191 .0787 .1228 11–20 .0651 .1104 .0629 .1053 .0656 .1105 .0668 .1069 40 5–10 .0666 .1102 .0716 .1195 .0696 .1134 .0740 .1198 11–20 .0617 .1063 .0631 .1097 .0699 .1129 .0662 .1078 50 5–10 .0704 .1142 .0682 .1152 .0685 .1128 .0722 .1181 11–20 .0612 .1093 .0677 .1115 .0582 .1036 .0694 .1135 100 5–10 .0634 .1102 .0662 .1148 .0661 .1153 .0654 .1162 11–20 .0597 .1074 .0605 .1086 .0610 .1062 .0588 .1062 200 5–10 .0575 .1115 .0590 .1082 .0592 .1084 .0595 .1110 11–20 .0596 .1088 .0560 .1061 .0577 .1025 .0633 .1092 300 5–10 .0588 .1081 .0553 .1024 .0641 .1138 .0595 .1096 11–20 .0549 .1067 .0527 .1011 .0587 .1071 .0590 .1058 400 5–10 .0544 .1001 .0572 .1075 .0565 .1013 .0622 .1100 11–20 .0546 .1045 .0569 .1059 .0590 .1023 .0546 .1046 500 5–10 .0574 .1080 .0584 .1081 .0610 .1089 .0565 .1046 11–20 .0545 .1029 .0516 .1004 .0573 .1085 .0590 .1076 Table 2 Observed proportion of rejections in 2,000 simulated data sets for each scenario, for functional data, under alternatives. The size of each sample has been randomly generated from a discrete uniform random law UD{a,a+1,...,b}, with (a,b)=(5,10), (11,20). The data were generated from a Wiener process where the mean function for population iis µi(t)=i kI(t>t0), 1 ≤i≤k,for several values of t0and different dispersion values σ. The nominal significance levels are α=0.05 and 0.10. Number of populations k t0Sample sizes limits ni Dispersion parameters of the Wiener processes σi=sσi∼U(1,s) s=1s=1.5s=2s=1.5s=2s=3 5% 10% 5% 10% 5% 10% 5% 10% 5% 10% 5% 10% 100 0.5 5–10 .9860 .9955 .6365 .7435 .3335 .4480 .8255 .9025 .5725 .6990 .2555 .3515 11–20 1.000 1.000 .9695 .9870 .7105 .8050 .9985 .9995 .9425 .9740 .5845 .7070 0.7 5–10 .8420 .9025 .3620 .4930 .1760 .2700 .4990 .6160 .3230 .4470 .1700 .2450 11–20 .9985 .9995 .7725 .8620 .4010 .5400 .9350 .9720 .7125 .8165 .3515 .4595 0.9 5–10 .2975 .4295 .1295 .2260 .0950 .1620 .1745 .2720 .1355 .2030 .0955 .1575 11–20 .6750 .7835 .2480 .3685 .1495 .2360 .3845 .5065 .2260 .3355 .1190 .1985 200 0.5 5–10 .9995 1.000 .8665 .9255 .5005 .6300 .9735 .9870 .7750 .8635 .4005 .5270 11–20 1.000 1.000 1.000 1.000 .9020 .9475 1.000 1.000 .9995 1.000 .7950 .8730 0.7 5–10 .9805 .9950 .5300 .6610 .2790 .3950 .7650 .8605 .4955 .6175 .2105 .3110 11–20 1.000 1.000 .9545 .9830 .6210 .7505 .9990 .9995 .9445 .9755 .4950 .6270 0.9 5–10 .4240 .5645 .1595 .2560 .1185 .2020 .2525 .3575 .1640 .2490 .0925 .1720 11–20 .9060 .9540 .3695 .5010 .1890 .2940 .5985 .7305 .3365 .4795 .1665 .2525 300 0.5 5–10 1.000 1.000 .9420 .9690 .6425 .7435 .9960 1.000 .9065 .9560 .4830 .6305 11–20 1.000 1.000 1.000 1.000 .9740 .9905 1.000 1.000 1.000 1.000 .9240 .9640 0.7 5–10 .9990 .9995 .6740 .7925 .3550 .4740 .8920 .9495 .6200 .7410 .2625 .3885 11–20 1.000 1.000 .9935 .9960 .7645 .8450 1.000 1.000 .9810 .9920 .6120 .7440 0.9 5–10 .5800 .7060 .2245 .3355 .1225 .2000 .3110 .4535 .1995 .3115 .1110 .1840 11–20 .9885 .9960 .4975 .6335 .2575 .3695 .7265 .8245 .4415 .5870 .1960 .2960 corresponds to a province. We have applied the proposed test to check the equality of the 52 mean functions and got a p-value <0.001 which leads us to reject the null hypothesis. The second data set consists of samples of Auslan (Australian Sign Language) signs. Auslan is the language used by the Australian Deaf and non-vocal communities. This data set is available in the Rpackage mfds (Górecki and Smaga [13]). In particular, 27 examples of each of 95 Auslan signs were captured from a native signer using high-quality position trackers and instrumented gloves and were collected over a period of nine weeks. The average length of each sign was approximately 57 frames. Each hand generated a total of 11 features: 3 for orientation (roll, pitch, yaw), 3 for position (x, y, z) and 5 for finger bends. Therefore, we have k=95 populations, 22 variables, and the sample sizes are n1=n2= ··· = n95 =27. For a more detailed description of the variables we refer to Kadous and Sammut [23]. In Fig. 2 the sample curves of the 2565 signs and the group sample means for the variable that measures the x-position of the left hand are represented. We have applied the test to check the equality of the 95 mean functions for each of the 22 variables. In all the cases the null hypothesis is rejected. 8 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 Fig. 1. All curves (left) and sample mean of each group (right) of the occupancy rate per month in Spain, from January 2011 to December 2019. There are k=52 populations (the 52 Spanish province), each sample curve is observed in 12 points (months) and the sample sizes are n1=n2= ··· = n52 =9. Fig. 2. All curves (left) and sample mean of each group (right) of the x-position of the left hand in the Auslan data set. There are k=95 populations (the different Auslan signs), 22 variables, and the sample sizes are n1=n2= ··· = n95 =27. 7. Comparison with existing tests In the Introduction we cited some tests for the comparison of the means of kindependent populations, that were designed for kfixed and increasing sample sizes. One may wonder if such tests are still valid in the setting considered in this paper. This section discusses such issue for the tests in Cuesta-Albertos and Febrero-Bande [8], Horváth and Rice [20], Cuevas et al. [9], Zhang et al. [31,32] and the approaches reviewed in Zhang [30]. The test in Cuesta-Albertos and Febrero-Bande [8] consists of randomly generating an element of L2(according to some adequate distribution), and then projecting the data on this element. The projected data has dimension 1, and so one can apply any technique designed for the univariate case. Notice that a random projection may lead to the rejection of the null hypothesis, while another random projection could induce the opposite conclusion. To avoid this inconvenience, these authors have proposed to take several random projections, calculate the p-value for each projection, and then apply some correction, as for example the procedure in Benjamini and Yekutieli [3], which controls the false discovery rate. In order to numerically investigate the validity of this test for large k, we partially repeated the experiment in Table 1. To apply the test in [8] we used the function fanova.RPm from the Rpackage fda.usc. For each case, we took 5, 10 and 30 random projections and adjusted the global p-value by using the procedure in [3]. Results based on 1000 simulated functions for each scenario are reported in Table 3. Looking at this table we see that the test becomes more liberal as kincreases, specially for smaller sample sizes. In the light of the outputs of this experiment, it can be concluded that the test in [8] is not valid in the setting of a large kand small ni. The test in Horváth and Rice [20] is also based on projections. Specifically, it considers the projections of the data on the eigenfunctions associated with the dlargest eigenvalues of an estimator of a certain covariance operator. These projections are also known as the scores. This way each function in the data set is transformed into a point of Rd. Then, these authors build a test statistic that compares the sample mean of the scores of each sample with an estimator of the 9 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 From Assumption 2, E{(∑1≤u=v≤ni⟨Xiu −µi,Xiv−µi⟩)4}≤n4 iM.(14) From (14),Assumptions 1 and 3and taking into account that k2var(Tk,Lin)≥∑k i=11 ni(ni−1) θi, 1 k4∑k i=11 n4 i(ni−1)4E{(∑1≤u=v≤ni⟨Xiu −µi,Xiv−µi⟩)4} var2(Tk,Lin)≤Mn4 max n4 min 1 k=o(1).(15) We have that E⎧ ⎨ ⎩(ni ∑ u=1⟨Xiu −µi, µi− ¯µ.⟩)4⎫ ⎬ ⎭=niE(⟨Xi−µi, µi− ¯µ.⟩4)+3ni(ni−1)E2(⟨Xi−µi, µi− ¯µ.⟩2) ≤4n2 iE(⟨Xi−µi, µi− ¯µ.⟩4).(16) From (16) and Assumption 4, 1 k4∑k i=11 n4 i E{(∑ni u=1⟨Xiu −µi, µi− ¯µ.⟩)4} var2(Tk,Lin)≤41 k 1 k∑k i=1E(⟨Xi−µi, µi− ¯µ.⟩4)/n2 i (1 k∑k i=1γi/ni)2=o(1).(17) Finally, (13),(15) and (17) imply (12), and hence the result is proven. □ Proof of Proposition 2.Let ε > 0, by Markov inequality, P(|ˆσ2 0k−σ2 0k|> σ2 0kε)≤E{(ˆσ2 0k−σ2 0k)2}/ε2σ2 0k.(18) Since E(ˆ θ2 i)=θ2 i, 1 ≤i≤k, we have that E{(ˆσ2 0k−σ2 0k)2} = (1/k2) k ∑ i=1 E{(ˆ θi−θi)2}/n2 i(ni−1)2.(19) From Assumption 2, routine calculations show that E{(ˆ θi−θi)2} ≤ M,∀i.(20) From (18)–(20) and Assumption 1, we get that P(|ˆσ2 0k−σ2 0k| σ2 0k > ε)≤M ε2 n4 max n4 min 1 k,(21) From Assumption 3, the right-hand side of (21) is o(1), which implies the result. □ Proof of the consistency ˆ θin the homoscedastic case. We have that |ˆ θ−θ| = ∫∫ |ˆ c(t,s)−c(t,s)||ˆ c(t,s)−c(t,s)|dtds. Therefore E(|ˆ θ−θ|)≤∫∫ E1/2[{ˆ c(t,s)−c(t,s)}2]E1/2[{ˆ c(t,s)+c(t,s)}2]dtds. Under the stated assumptions, routine calculations show that E[{ˆ c(t,s)−c(t,s)}2] ≤ C/kand E[{ˆ c(t,s)+c(t,s)}2] ≤ C, for all 0 ≤t,s≤1, for a positive constant C, which implies that E(|ˆ θ−θ|)→0, and hence the consistency of ˆ θas an estimator of θ.□ Proof of Proposition 3.We will show it by reduction to absurdity. Suppose that 1 σ2 0k 1 k∑k i=1 γi ni ↛0 which is equivalent to ∃ε1>0 such that ε1<1 σ2 0k 1 k k ∑ i=1 γi ni for an infinite number of values of k. (22) 16 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 By applying the Cauchy–Schwarz inequality and Assumption 2, one gets 1 σ2 0k 1 k k ∑ i=1 γi ni≤MDk σ2 0k .(23) From (22) and (23), 0< ε2=ε1 M≤Dk σ2 0k ,for an infinite number of values of k. (24) From Assumption 1 it follows that σ2 0k≥τ/n2 max (25) Let ε3>0 be arbitrary but fixed. We are assuming that √kDk σk→δ∈ [0,∞), which implies that ∃k0=k0(ε3) such that √kDk σk≤δ+ε3,∀k≥k0.(26) From (23),(25) and (26), M(k n2 max )1/2Dk/σ2 0k √1+Dk/σ2 0k ≤δ+ε3,∀k≥k0.(27) Taking into account that the function f(x)=x/√1+xis increasing ∀x≥0, from (24) it follows that Dk/σ2 0k √1+Dk/σ2 0k ≥ε2 √1+ε2 ,for an infinite number of values of k. (28) From (27) and (28) it follows that M(k/n2 max)1/2≤δ+ε3, for an infinite number of values of k, which contradicts (9), implying that 1 σ2 0k 1 k∑k i=1 γi ni→0. □ Proof of Theorems 3 and 4.Their proofs are parallel to that of Lemmas 1 and 2. To save space we omit it. □ Acknowledgments The authors thank two anonymous referees for their valuable time and careful comments, which improved the presentation of this paper. The authors have been partially supported by grant MTM2017-89422-P (Spanish Ministry of Economy, Industry and Competitiveness, the State Agency of Investigation, the European Regional Development Fund). MD Jiménez-Gamero was also partially supported by grant P18-FR-2369 (Junta de Andalucía, Spain). Appendix A. Supplementary data Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jmva.2021.104778. References [1] M.G. Akritas, S. Arnold, Asymptotics for analysis of variance when the number of levels is large, J. Amer. Statist. Assoc. 95 (2000) 212–226. [2] M.G. Akritas, N. Papadatos, Heteroscedastic one-way ANOVA and lack-of-fit tests, J. Amer. Statist. Assoc. 99 (2004) 368–382. [3] Y. Benjamini, D. Yekutieli, The control of the false discovery rate in multiple testing under dependency, Ann. Statist. 29 (2001) 1165–1188. [4] M. Benko, W. Härdle, A. Kneip, Common functional principal components, Ann. Statist. 37 (2009) 1–34. [5] D.D. Boos, C. Brownie, ANOVA and rank tests when the number of treatments is large, Statist. Probab. Lett. 23 (1995) 183–191. [6] M. Cousido-Rocha, J. de Uña Álvarez, J.D. Hart, Testing equality of a large number of densities under mixing conditions, Test 28 (2019) 1203–1228. [7] J.A. Cuesta-Albertos, E. del Barrio, R. Fraiman, C. Matrán, The random projection method in goodness of fit for functional data, Comput. Statist. Data Anal. 51 (2007) 4814–4831. [8] J.A. Cuesta-Albertos, M. Febrero-Bande, A simple multiway ANOVA for functional data, Test 19 (2010) 537–557. [9] A. Cuevas, M. Febrero, R. Fraiman, An ANOVA test for functional data, Comput. Statist. Data Anal. 47 (2004) 111–122. [10] B. Efron, C. Stein, The jackknife estimate of variance, Ann. Statist. 9 (1981) 586–596. [11] A. Ghiglietti, F. Ieva, A.M. Paganoni, Statistical inference for stochastic processes: two-sample hypothesis tests, J. Statist. Plann. Inference 180 (2017) 49–68. [12] T. Górecki, S. Hörmann, L. Horváth, P. Kokoszka, Testing normality of functional time series, J. Time Anal. 39 (2018) 471–487. [13] T. Górecki, L. Smaga, mfds: Multivariate Functional Data Sets. R package version 0.1.0, 2017, URL https://github.com/Halmaris/mfds. [14] T. Górecki, L. Smaga, fdANOVA: an R software package for analysis of variance for univariate and multivariate functional data, Comput. Statist. 34 (2019) 571–597. [15] A.K. Gupta, S.W. Harrar, Y. Fujikoshi, MANOVA for large hypothesis degrees of freedom under non-normality, Test 17 (2008) 120–137. 17 M.D. Jiménez-Gamero and A.M. Franco-Pereira Journal of Multivariate Analysis 185 (2021) 104778 [16] S.W. Harrar, A.C. Bathke, A modified two-factor multivariate analysis of variance: asymptotics and small sample approximations, Ann. Inst. Statist. Math. 64 (2012) 135–165. [17] S.W. Harrar, A.K. Gupta, Asymptotic expansion for the null distribution of the F-statistic in one-way ANOVA under non-normality, Ann. Inst. Statist. Math. 59 (2007) 531–556. [18] N. Henze, M.D. Jiménez-Gamero, A test for Gaussianity in Hilbert spaces via the empirical characteristic functional, Scand. J. Stat. (2020) http://dx.doi.org/10.1111/sjos.12470. [19] L. Horváth, P. Kokoszka, Inference for Functional Data with Applications, Springer, New York, 2012. [20] L. Horváth, G. Rice, An introduction to functional data analysis and a principal component approach for testing the equality of mean curves, Rev. Mat. Complut. 28 (2015) 505–548. [21] Q. Jiang, M. Hušková, S.G. Meintanis, L. Zhu, Asymptotics, finite-sample comparisons and applications for two-sample tests with functional data, J. Multivariate Anal. 170 (2019) 202–220. [22] M.D. Jiménez-Gamero, M. Cousido-Rocha, M.V. Alba-Fernández, F. Jiménez-Jiménez, Testing the equality of a large number of populations, TEST (2021) http://dx.doi.org/10.1007/s11749-021-00769-9. [23] M.W. Kadous, C. Sammut, Classification of multivariate time series and structured data using constructive induction, Mach. Learn. 58 (2005) 179–216. [24] E. Paparoditis, T. Sapatinas, Bootstrap-based testing of equality of mean functions or equality of covariance operators for functional data, Biometrika 103 (2016) 727–733. [25] J. Park, D. Park, Testing the equality of a large number of normal population means, Comput. Statist. Data Anal. 56 (2012) 1131–1149. [26] R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2020, URL https://www.R-project.org/. [27] R.J. Serfling, Approximation Theorems of Mathematical Statistics, Wiley, New York, 1980. [28] A. Yuan, H.B. Fang, H. Li, C.O. Wu, M.T. Tan, Hypothesis testing for multiple mean and correlation curves with functional data, Statist. Sinica 30 (2020) 1095–1116. [29] D. Zhan, J.D. Hart, Testing equality of a large number of densities, Biometrika 101 (2014) 449–464. [30] J.T. Zhang, Analysis of Variance for Functional Data, Chapman & Hall, London, 2013. [31] J.T. Zhang, M.Y. Cheng, H.T. Wu, B. Zhou, A new test for functional one-way ANOVA with applications to ischemic heart screening, Comput. Statist. Data Anal. 132 (2019) 3–17. [32] J.T. Zhang, X. Liang, One-way ANOVA for functional data via globalizing the pointwise F-test, Scand. J. Stat. 41 (2014) 51–71. [33] C. Zhang, H. Peng, J.T. Zhang, Two samples tests for functional data, Commun. Stat. - Theory Methods 39 (2010) 559–578. 18