Two-sample Tests for Functional Data Using Characteristic Functions

In this paper, we consider the two-sample problem for univariate and multivariate functional data. To solve this problem, we use tool of characteristic function and a basis function representation of functional data. We construct test statistics for conformity of distributions based on a weighted distance between characteristic functions of random vectors obtained in basis representation. Different weight functions result in different test statistics, whose distributions are approximated by permutation method. Testing procedures are implemented in the R program and the code is available. Simulation study shows good finite sample properties of proposed methods, while real data example illustrates the application of them.


Introduction
Over the last three decades, real-time measurement instruments and data storage computer resources have significantly developed. This gives possibility to observe and save large amount of data as results of random experiments. In particular, very dense continuous-time monitoring clinical diagnostics or stock market information or meteorological information are common nowadays. Such data can be seen as functions of time or space. In Functional Data Analysis (FDA), they are in fact treated in such a way, i.e., as samples of random functions. Functional data can be univariate, when we observe one process (e.g., temperature), or multivariate, when we consider multiple processes simultaneously (e.g., temperature and precipitation). They can also be dense or sparse. FDA deals with the statistical description and modelization of samples of this type. Most of statistical methods of standard statistic have their counterparts for functional data. For details, we refer to the following monographs and the references given therein: Ferraty and Vieu (2006); Horváth and Kokoszka (2012); Silverman (2002, 2005); Zhang (2013). Many FDA methods are implemented in the R program (R Core Team 2019). Perhaps the most well known and the biggest packages for FDA are the fda (Ramsay, Wickham, Graves, and Hooker 2018) and fda.usc (Febrero-Bande and Oviedo de la Fuente 2012) packages. For review of R packages for FDA, see the supplementary materials to the paper Górecki and Smaga (2019).
The two-sample test for conformity of distributions is one of the statistical methods that have been generalized for functional data. However, it is rarely considered in the literature, and it is not so popular as the analysis of variance problem for functional data (see, for example, Smaga 2017, 2019;Zhang 2013). The main two papers considering the two-sample problem for functional data are Hall and Van Keilegom (2007) and Pomann, Staicu, and Ghosh (2016). Nevertheless, we propose other approach to solve the problem than those in these papers, so we do not consider them in details. The present paper considers methods based on characteristic functions of random processes and their basis function representation. This representation can be seen as dimension reduction method, which represents functional data as random vectors of much less dimension than discrete functional data. The two-sample tests for multivariate data using characteristic functions have been already considered in the literature. Perhaps the most famous one is based on the so-called energy statistic proposed by Székely and Rizzo (2013). This test statistic is in fact an estimator of weighted distance between characteristic functions of two random vectors with certain weight function. Other weight functions can be also considered. Recently, Chen, Meintanis, and Zhu (2019) investigated this problem and proposed to use a density of some random distribution as the weight function. Their method has similar theoretical properties to the energy statistic. Moreover, it gives the opportunity to use different weight functions, i.e., different densities. For instance, Chen et al. (2019) used densities of spherical stable distributions, which greatly simplified test statistic. In this paper, we extend these methods for univariate and multivariate functional data. First, we express the two-sample problem for these data in terms of characteristic functions of random processes. Then using the basis function representation, we reduce this problem to the twosample problem for random vectors. Finally, we propose permutation tests constructed based on results of Székely and Rizzo (2013) and Chen et al. (2019). Simulation study and real data example show that the new tests maintain the type I error level very well and have satisfactory power. However, there is not one test, which performs best. Nevertheless, the tests using density weight function are usually more powerful than test using energy statistic.
The remainder of this paper is as follows. In Section 2, we formulate the two-sample problem for univariate and multivariate functional data and present its solution using characteristic functions, basis function representation and methods for two-sample problem for random vectors. Implementation of the new tests is also discussed there. Section 3 includes the Monte Carlo study on the finite sample behavior of proposed methods in terms of size control and power. In Section 4, a real data illustration is given. Finally, the conclusions are presented in Section 5. The supplementary materials to this paper contain the code of R programming language (R Core Team 2019) for numerical experiments.

Two-sample problem for functional data
In this section, we present the two-sample problem for functional data and its alternative form in the framework of characteristic functions. Next using a basis function representation of functional data, we consider a weighted distance between characteristic functions being a base to construct test statistics. Different weight functions are also discussed. The permutation tests for two-sample problem for functional data are constructed based on estimators of weighted distances. Finally, we discuss the implementation of proposed testing procedures.

Hypotheses and base for test statistic
, where a, b ∈ R and a < b, be the Hilbert space of p-dimensional vectors of square integrable functions defined on the interval [a, b] and endowed with the following inner product: a, b]. We assume that X 1 and X 2 are two independent random processes belonging to the Hilbert space L p 2 [a, b], and we wish to test the following null hypothesis: where d = stands for equality in distribution. To solve this problem, we would like to use the fact that the null hypothesis H 0 holds if and only if the characteristic functions of the distributions of processes X 1 and X 2 are equal. By (Bosq 2000, p. 37), the characteristic functions of X 1 and X 2 are as follows ϕ X j (u) = E(exp(i u, X j p )) for j = 1, 2, u ∈ L p 2 [a, b] and i 2 = −1. Note that we suppose that for all u ∈ L p 2 [a, b] the integrals u, X j p converge for almost all realizations of X j , j = 1, 2.
To effectively deal in practice with the characteristic function given above, we use the basis function representation of random processes X 1 and X 2 . Let {φ kl } ∞ l=1 be basis in L 1 2 [a, b], k = 1, . . . , p, and X j = (X j1 , . . . , X jp ) for j = 1, 2. Then where α jkl are random variables with finite variance, j = 1, 2, k = 1, . . . , p. However, it is not possible to use such representation in practice. Fortunately, Ramsay and Silverman (2005) noticed that only a number of the first coefficients in this representation is usually the largest and the most practically important. For this reason, we use the finite representation as follows: assuming that the processes X jk can be represented by a finite (appropriately large) number B k of basis functions. This representation will be called the basis function representation of random process. In the following, we will use the matrix notation of basis function representation given in Equation (1). Namely where α j = (α j11 , . . . , α j1B 1 , . . . , α jp1 , . . . , α jpBp ) ∈ R B are random vectors j = 1, 2, Some practical remarks about the basis function representation are as follows (see, for example, Krzyśko and Waszak 2013). The coefficients α jkl are estimated for each coordinate separately using the least squares method or the regularized maximum likelihood method (Matsui, Araki, and Konishi 2008) or the roughness penalty approach (Ramsay and Silverman 2005). The values of B k may be selected depending on the aim of the research, for instance, for the best fit, the Bayesian information criterion should perhaps be used (Shmueli 2010). Note that these values determine the degree of smoothness. Namely, small value of B k causes more smoothness. In practice, it is usually not crucial, which basis is used. However, the Fourier basis is usually suggested for periodic or nearly periodic data, while the B-spline basis is typically used for non-periodic locally smooth data (Horváth and Kokoszka 2012).
Let us return to the characteristic functions of random processes X 1 and X 2 . Assuming that the function where γ = J Φ β, which means that the characteristic functions of the random processes X 1 and X 2 are the characteristic functions of the random vectors α 1 and α 2 respectively. Thus, to test the null hypothesis H 0 , we propose to use methods based on the following distance between characteristic functions for random vectors: where |z| is the modulus of z ∈ C, and w > 0 (almost everywhere) is a weight function. The weight function can be chosen in many ways. In the next section, we discuss two approaches considered in the literature.

Weight functions
Well known statistic of the form given in Equation (2) is the two-sample energy statistic of Székely and Rizzo (2013), which uses the weight function w = e, where and · is the standard Euclidean norm in the space R B . Then Székely and Rizzo (2013) proved that D e = 0 if and only if α 1 d = α 2 , and it is otherwise strictly positive. This implies that testing the null hypothesis H 0 is equivalent to testing H De 0 : D e = 0. The test for this hypothesis is constructed in the next section.
On the other hand, Chen et al. (2019) proposed to use a density of some random variable as the weight function. Let f be a density (positive with probability 1) and let ϕ be a corresponding characteristic function of a B × 1 random vector. Then for w = f , we have where Re(z) denotes the real part of z ∈ C. Furthermore, Chen et al. (2019) showed that D f = 0 if and only if α 1 d = α 2 , and it is otherwise strictly positive. Thus, similarly to the energy statistic, we can test H 0 by testing H We can use different choices of the density f . Chen et al. (2019) considered the densities f α of spherical stable distributions with exponents α ∈ (0, 2], which greatly simplify D f . The characteristic functions of these distributions are as follows: For α = 1 and α = 2, we have the multivariate standard Cauchy and normal distributions respectively. Taking f = f α , we obtain which we will use.

Estimation and permutation test
In practice, we first have to estimate D e and D fα based on the data. Let X 11 , . . . , X 1n 1 and X 21 , . . . , X 2n 2 be independent samples for random processes X 1 and X 2 respectively. We represent these observations in basis as presented in Section 2.1, i.e., X ji (t) = Φ(t)α ji for j = 1, 2, i = 1, . . . , n j . Then, the natural estimators of D e and D fα are as followŝ : D fα = 0, we propose permutation tests based on test statisticsD e andD fα respectively. We reject the null hypotheses for large values of these test statistics. The permutation test based on test statisticD ∈ {D e ,D fα } has the following steps: 1. ComputeD for original data α ji , j = 1, 2, i = 1, . . . , n j .
2. Create a permutation sample from the given data in the following way: From all observations α ji , j = 1, 2, i = 1, . . . , n j , select randomly without replacement n 1 observations for the first new sample, and the remainder of the observations creates the second sample.
5. The final p-value of the permutation test is defined by P −1 P m=1 I(D m >D), where I(A) denotes the indicator function on the set A (takes value 1 if A is true and 0 otherwise).

Implementation
The permutation two-sample tests for functional data proposed above, the simulation experiments of Section 3 as well as real data example of Section 4 were implemented in the R program (R Core Team 2019). The code is available in the supplementary materials. To obtain a basis function representation of functional data, the following functions from the fda package (Ramsay et al. 2018) were used: create.fourier.basis(), create.bspline.basis() and smooth.basis(). The implementation of permutation test based on two-sample energy statis-ticD e in the function eqdist.etest() available in the energy package (Rizzo and Székely 2019) was applied. The permutation tests based on test statisticsD fα were implemented by the authors using the Rcpp package (Eddelbuettel and François 2011;Eddelbuettel 2013;Eddelbuettel and Balamuta 2017), which significantly improves the execution time.

Simulation studies
In this section, we describe the simulation experiments and their results, in which we study the finite sample properties of the tests proposed in Section 2. Namely, we consider the empirical size and power of the permutation tests based on test statisticsD e ,D f 0.1 ,D f 0.5 ,D f 1 ,D f 1.5 andD f 2 .

Experimental setup
To calculate the values of test statistics, we used the Fourier and B-spline bases with B k = 5 for k = 1, . . . , p. The coefficients of basis function representation were estimated by the least squares method. We used 1, 000 simulation and permutation samples to estimate the empirical sizes and powers and p-values respectively. For simplicity, the significance level was set to 5%.
We set [a, b] = [0, 1]. The discrete functional data for the above processes were generated in an equally spaced grid of m design time points in the interval [0, 1]. Namely, the Wiener and Ornstein-Uhlenbeck processes were observed in points t 1 = 0, . . . , t m+1 = 1. However, the finale observations used were those for m points t 2 , . . . , t m+1 , since we removed the first zero value of these processes for t 1 . On the other hand, the Brownian bridge was observed in points t 1 = 0, . . . , t m+2 = 1, while the finale observations used were those for m points t 2 , . . . , t m+1 , i.e., we removed the first and the last zero value of Brownian bridge for t 1 and t m+2 . We set m = 25, 50 for p = 1 and m = 25 for p = 3. For p = 3, additionally to the above independent case (the coordinates of both processes X 1 and X 2 are generated independently), we consider the dependent case described as follows. To all generated observations of processes X 1 and X 2 , we add a kind of random errors e ji , j = 1, 2, i = 1, . . . , n j , whose coordinates are dependent random processes. Let N ji ∼ N 9 (0 9 , Σ) be independent random vectors, where 0 a is the a × 1 vector of zeros, Σ = σ (1 − ρ)I 9 + ρ1 9 1 9 with σ = 0.2, ρ = 0.1, I a is the a × a identity matrix, and 1 a is the a × 1 vector of ones. Then e ji (t) = Φ(t)N ji , where t ∈ [0, 1], the 3 × 9 matrix Φ is as in Section 2 and contains the Fourier basis functions with B k = 3 for k = 1, 2, 3.

Simulation results
Let us discuss the simulation results, which are presented in Tables 1-2 and Table 3 for univariate and multivariate cases respectively.
First of all, we can observe that the empirical sizes (δ = 0) of all permutation tests are very close to the nominal level of significance 5% in all scenarios considered. Thus all tests maintain the type I error level very well. This follows from that the distributions of both processes X 1 and X 2 under the null hypothesis are the same, and hence the permutation tests (using all permutations) should be exact, i.e., the type I error level is equal to the significance level. Now we consider the power (δ > 0). The empirical powers of all tests increase with the increase of the number of observations (compare, for example, the results in Tables 1 and 2). They also increase, when δ increases, i.e., "we are moving away from the null hypothesis". The empirical power of all tests seems to not depend on the number of design time points, at which the functional data are observed (Tables 1-2, m = 25, 50). We also note that under dependent scenario, the empirical power of all tests is less than this for independent one, but this is not surprising as in dependent case the observations of original processes (Wiener, etc.) were contaminated by random noise.
In univariate case (Tables 1-2), theD fα tests are much more powerful than theD e test. We Table 1: Empirical sizes (δ = 0) and powers (δ > 0) (as percentages) of all tests obtained for p = 1 and n 1 = n 2 = 15. Column B denotes basis used (F -Fourier, B -B-spline). observe very similar behavior in the multivariate case (Table 3), but for some cases with B-spline basis in independent scenario, the empirical power of theD e test is greater than this of the worst (in term of power) of theD fα tests. The power comparison for theD fα tests with respect to α is not easy. There is not oneD fα test, which has the greatest power in all cases. In one-dimensional case with Fourier basis, the empirical powers of theD fα tests are usually similar (under Wiener and Ornstein-Uhlenbeck processes) or decrease with increase of α (under Brownian bridge). On the other hand, in multi-dimensional and independent case with Fourier basis, the empirical powers of theD fα tests increase with increase of α. When using B-spline basis or in dependent multivariate case, the situation changes. In these cases, theD fα tests with α = 1, 1.5, 0.5 usually overcome the remaining tests.

Real data example
In this section, we illustrate the practical use of the permutation tests proposed in Section 2. For this purpose, we consider the Berkeley Growth Study data Silverman 2002, 2005;Tuddenham and Snyder 1954). The data are available in the data set growth contained in the R package fda (Ramsay et al. 2018). This data set contains the heights in centimeters of 39 boys and 54 girls measured at a set of 31 ages from age 1 to 18. The ages are not equally spaced; there are four measurements while the child is one year old, annual measurements from two to eight years, followed by heights measured biannually. Thus the heights can be treated as values of some random function of age and the data set as sample of 93 univariate (p = 1) functional observations measured at 31 design time points. The trajectories of the growth curves for boys and girls separately are presented in Figure 1.
The growth data are naturally divided into two groups, namely boys and girls. Thus we have two samples of functional data of sizes n 1 = 39 and n 2 = 54 respectively. For illustrative purposes, we would like to test the equality of distributions of growth of boys and girls, which seems to be false, since the height of boys is usually greater than the height of girls. For this problem, we applied the permutation tests based on test statisticsD e ,D f 0.1 ,D f 0.5 ,D f 1 , D f 1.5 andD f 2 using 1, 000 permutation samples. Since these data seem to be non-periodic, we use the B-spline basis only with B 1 = 5, 6, . . . , 15. The p-values of all testing procedures are presented in Table 4.
We observe that theD e ,D f 0.1 ,D f 0.5 andD f 1 tests have p-values equal to zero or close to zero, and hence they reject the null hypothesis. For the remaining tests, the situation changes.
Their p-values are usually much larger. For small (B 1 = 5, 6) and larger (B 1 = 15) numbers of basis functions used in basis representation, these testing procedures do not reject the null hypothesis, while for the other values of B 1 (i.e., B 1 = 7, . . . , 14), they reject the null hypothesis. To indicate which decision for these tests should be made, we note that optimal value of B 1 according to the procedure in Krzyśko and Waszak (2013) using the Bayesian information criterion (BIC) is equal to 11. Thus perhaps we should use B 1 around 11, which means that the tests reject the null hypothesis. Therefore the new tests confirm the fact that the boys are usually higher than girls.  Table 4: P-values of all tests using B-spline basis for testing equality of distributions of heights of boys and girls contained in the data set growth. The bold row concerns B 1 = 11, which is optimal number of basis functions according to BIC.

Conclusions
We have proposed the new permutation tests for the two-sample problem for functional data. Our approach covers both univariate and multivariate cases. Using characteristic function representation of this two-sample problem and basis function representation of functional data, we have reduced this problem to the multivariate two-sample problem. Then we have applied the tests for the last problem, which are based on weighted distances between characteristic functions as test statistics. Simulation study and real data application have been conducted for determining the finite sample properties of new tests. They have showed that the permutation tests control the type I error level very well and have satisfactory power. Moreover, in term of power, the tests using density weight function are better than test using energy statistic in most cases. The proposed method may be extended, for instance, by using other weight functions to construct the test statistic.