Asymptotic and bootstrap tests for the dimension of the non-Gaussian subspace

TL;DR

This paper introduces asymptotic and bootstrap statistical tests to accurately determine the dimension of the non-Gaussian subspace in data, enhancing the effectiveness of non-Gaussian component analysis.

Contribution

It develops new asymptotic and bootstrap testing procedures for the subspace dimension in non-Gaussian component analysis using FOBI.

Findings

Tests perform well in simulations

Bootstrap method improves accuracy in finite samples

Provides reliable dimension estimation in practice

Abstract

Dimension reduction is often a preliminary step in the analysis of large data sets. The so-called non-Gaussian component analysis searches for a projection onto the non-Gaussian part of the data, and it is then important to know the correct dimension of the non-Gaussian signal subspace. In this paper we develop asymptotic as well as bootstrap tests for the dimension based on the popular fourth order blind identification (FOBI) method.

Figures

Tables

Table 1. TABLE I: Model M1: Rejection rates in 1000 repetitions for asymptotic and bootstrap tests of H 0 , 1 subscript 𝐻 0 1 H_{0,1} , H 0 , 2 subscript 𝐻 0 2 H_{0,2} (true) and H 0 , 3 subscript 𝐻 0 3 H_{0,3} with α = 0.05 𝛼 0.05 \alpha=0.05 .

	$H_{0,1}$		$H_{0,2}$		$H_{0,3}$
$n$	Asy	Boot	Asy	Boot	Asy	Boot
500	0.040	0.128	0.003	0.031	0.001	0.012
1000	0.149	0.252	0.004	0.058	0.000	0.016
2000	0.403	0.466	0.027	0.068	0.000	0.024
5000	0.950	0.950	0.040	0.059	0.001	0.011
10000	0.999	0.999	0.053	0.058	0.002	0.010

Table 2. TABLE II: Model M2: Rejection rates in 1000 repetitions for asymptotic and bootstrap tests of H 0 , 3 subscript 𝐻 0 3 H_{0,3} , H 0 , 4 subscript 𝐻 0 4 H_{0,4} (true) and H 0 , 5 subscript 𝐻 0 5 H_{0,5} with α = 0.05 𝛼 0.05 \alpha=0.05 .

	$H_{0,3}$		$H_{0,4}$		$H_{0,5}$
$n$	Asy	Boot	Asy	Boot	Asy	Boot
500	0.577	0.328	0.050	0.065	0.000	0.008
1000	0.997	0.964	0.073	0.056	0.002	0.008
2000	1.000	1.000	0.057	0.049	0.002	0.016
5000	1.000	1.000	0.068	0.064	0.004	0.018
10000	1.000	1.000	0.050	0.046	0.050	0.046

Table 3. TABLE III: Model M3: Rejection rates in 1000 repetitions for asymptotic and bootstrap tests of H 0 , 2 subscript 𝐻 0 2 H_{0,2} , H 0 , 3 subscript 𝐻 0 3 H_{0,3} (true) and H 0 , 4 subscript 𝐻 0 4 H_{0,4} with α = 0.05 𝛼 0.05 \alpha=0.05 .

	$H_{0,2}$		$H_{0,3}$		$H_{0,4}$
$n$	Asy	Boot	Asy	Boot	Asy	Boot
500	0.715	0.796	0.024	0.051	0.001	0.015
1000	0.993	0.995	0.020	0.034	0.000	0.013
2000	1.000	1.000	0.036	0.044	0.006	0.016
5000	1.000	1.000	0.042	0.041	0.003	0.011
10000	1.000	1.000	0.043	0.042	0.005	0.012

Equations

$$\mathbf{x}=\mathbf{A}\mathbf{z}+\mathbf{b}$$

$$\mathbf{x}={\mathbf{A}}_1{\mathbf{z}}_1+{\mathbf{A}}_2{\mathbf{z}}_2+\mathbf{b}$$

$${\mathbf{S}}_1=E\left[(\mathbf{x}-E(\mathbf{x}))(\mathbf{x}-E(\mathbf{x}){)}^{\prime }\right]$$

$${\mathbf{S}}_2=E\left[r^2(\mathbf{x}-E(\mathbf{x}))(\mathbf{x}-E(\mathbf{x}){)}^{\prime }\right],$$

$$\mathbf{W}{\mathbf{S}}_1{\mathbf{W}}^{\prime }={\mathbf{I}}_p\text{mbox}and\mathbf{W}{\mathbf{S}}_2{\mathbf{W}}^{\prime }=\mathbf{D}$$

$$\mathbf{R}={\mathbf{S}}_1^{-1/2}{\mathbf{S}}_2{\mathbf{S}}_1^{-1/2}.$$

$$m_1(\mathbf{R}):=tr(\mathbf{R})/p\text{mbox}and{m}_2(\mathbf{R}):=tr({\mathbf{R}}^2)/p,$$

$$H_{0,k}:\ \ \mbox{exactly $p-k$ eigenvalues of $\boldsymbol{R}$ equal $p+2$ }$$

$${\hat{\mathbf{S}}}_1=\frac{1}{n}\sum _{i=1}^n({\mathbf{x}}_i-\bar{\mathbf{x}})({\mathbf{x}}_i-\bar{\mathbf{x}}{)}^{\prime },$$

$${\hat{\mathbf{S}}}_2=\frac{1}{n}\sum _{i=1}^n({\mathbf{x}}_i-\bar{\mathbf{x}})({\mathbf{x}}_i-\bar{\mathbf{x}}{)}^{\prime }{\hat{\mathbf{S}}}_1^{-1}({\mathbf{x}}_i-\bar{\mathbf{x}})({\mathbf{x}}_i-\bar{\mathbf{x}}{)}^{\prime }$$

$$T_k=\underset{\mathbf{U}\in {\mathcal{O}}^{p\times (p-k)}}{\min }{m}_2\left({\mathbf{U}}^{\prime }(\hat{\mathbf{R}}-(p+2){\mathbf{I}}_p)\mathbf{U}\right).$$

$${\hat{\mathbf{U}}}_k=\underset{\mathbf{U}\in {\mathcal{O}}^{p\times (p-k)}}{\mathrm{argmin}}{m}_2\left({\mathbf{U}}^{\prime }(\hat{\mathbf{R}}-(p+2){\mathbf{I}}_p)\mathbf{U}\right)$$

$$T_k^{\ast }=m_2\left({\mathbf{V}}_k^{\prime }(\hat{\mathbf{R}}-(p+2){\mathbf{I}}_p){\mathbf{V}}_k\right)$$

$$C_k\sim 2{\sigma }_1{\chi }_{(p-k-1)(p-k+2)/2}^2+\left(2{\sigma }_1+{\sigma }_2(p-k)\right){\chi }_1^2$$

$$n(p-k)T_k\ge c_{k,\alpha }$$

$$\mathbb{P}(n(p-k)T_{k}\geq c_{k,n})\ \rightarrow\ \left\{\begin{array}[]{ll}1,&\hbox{if $k<q$ ;}\\ 0,&\hbox{if $k\geq q$.}\end{array}\right.$$

$$\hat{q}=\min \{k:n(p-k)T_k<c_{k,n}\}{\to }_{P}q.$$

$${\mathbf{x}}_i^{\ast }={\hat{\mathbf{P}}}_k{\tilde{\mathbf{x}}}_i+{\hat{\mathbf{S}}}_1^{1/2}{\hat{\mathbf{U}}}_k{\mathbf{o}}_i,i=1,...,n,$$

$$\hat{p}=\frac{\#(T_i^{\ast }\ge T)+1}{M+1}.$$

$$\frac{n(p-k)T_{k,1}}{2{\hat{\sigma }}_1}\text{mbox}or\frac{n(p-k)T_{k,2}}{2{\hat{\sigma }}_1+4(p-k)}$$

$$T_k=\underset{\mathbf{U}\in {\mathcal{O}}^{p\times (p-k)}}{\min }{s}^2\left({\mathbf{U}}^{\prime }\hat{\mathbf{R}}\mathbf{U}\right).$$

$$\mathbf{x}=\mathbf{A}\mathbf{z}+\mathbf{b}$$

$${\tilde{\mathbf{S}}}_1=\frac{1}{n}\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }\text{mbox}and{\tilde{\mathbf{S}}}_2=\frac{1}{n}\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }{\mathbf{x}}_i{\mathbf{x}}_i^{\prime }$$

$$\hat{\mathbf{R}}={\hat{\mathbf{S}}}_1^{-1/2}{\hat{\mathbf{S}}}_2{\hat{\mathbf{S}}}_1^{-1/2}$$

$$\widehat{\boldsymbol{R}}=\left(\begin{array}[]{cc}\widehat{\boldsymbol{R}}_{11}&\widehat{\boldsymbol{R}}_{12}\\ \widehat{\boldsymbol{R}}_{21}&\widehat{\boldsymbol{R}}_{22}\\ \end{array}\right).$$

$$\sqrt{n}(\hat{\mathbf{R}}-(p+2){\mathbf{I}}_p{)}_{22}$$

$$\tilde{\mathbf{R}}$$

$${\mathbf{U}}^{\prime }{\mathbf{V}}_k^{\prime }(\tilde{\mathbf{R}}-(p+2){\mathbf{I}}_p){\mathbf{V}}_k\mathbf{U}\sim {\mathbf{V}}_k^{\prime }(\tilde{\mathbf{R}}-(p+2){\mathbf{I}}_p){\mathbf{V}}_k$$

$$\sqrt{n}({\mathbf{V}}_k^{\prime }(\hat{\mathbf{R}}-(p+2){\mathbf{I}}_p){\mathbf{V}}_k){\to }_d{\mathbf{N}}_k$$

$${\sigma }_1\left({\mathbf{I}}_{(p-k{)}^2}+{\mathbf{K}}_{p-k,p-k}\right)+{\sigma }_{2}vec({\mathbf{I}}_{p-k})vec({\mathbf{I}}_{p-k}{)}^{\prime },$$

Full text

I Introduction

Throughout the paper we assume the Non-Gaussian Component Analysis (NGCA) model, that is, $\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$ is a random sample from the distribution of

[TABLE]

where $E(\boldsymbol{z})=\boldsymbol{0}$ and $Cov(\boldsymbol{z})=\boldsymbol{I}_{p}$, $\boldsymbol{A}\in\mathbb{R}^{p\times p}$ is nonsingular, $\boldsymbol{b}\in\mathbb{R}^{p}$, and $\boldsymbol{z}=(\boldsymbol{z}_{1}^{\prime},\boldsymbol{z}_{2}^{\prime})^{\prime}$, $\boldsymbol{z}_{1}\in\mathbb{R}^{q}$ and $\boldsymbol{z}_{2}\in\mathbb{R}^{p-q}$ are independent random vectors, $\boldsymbol{z}_{1}$ is non-Gaussian and $\boldsymbol{z}_{2}$ is Gaussian. For non-Gaussian $\boldsymbol{z}_{1}$, there is no $\boldsymbol{a}\in\mathbb{R}^{q}$ such that $\boldsymbol{a}^{\prime}\boldsymbol{z}_{1}$ has a normal distribution. For a Gaussian $\boldsymbol{z}_{2}$, $\boldsymbol{a}^{\prime}\boldsymbol{z}_{2}$ has a normal distribution for all $\boldsymbol{a}\in\mathbb{R}^{p-q}$. The idea then is, based on the observations $\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$, to make inference on the unknown $q$, $0\leq q\leq p$, and estimate the non-Gaussian signal and Gaussian noise subspaces determined by $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$, respectively. In the literature it is often preassumed that the dimension $q$ is known. In this short note we develop tests for $q$ that are based on the matrix of fourth moments used in fourth order blind identification (FOBI).

The model can also be written as

[TABLE]

where now $\boldsymbol{A}_{1}\in\mathbb{R}^{p\times q}$ and $\boldsymbol{A}_{2}\in\mathbb{R}^{p\times(p-q)}$ have ranks $q$ and $p-q$, respectively. The independent random vectors $\boldsymbol{A}_{1}\boldsymbol{z}_{1}$ and $\boldsymbol{A}_{2}\boldsymbol{z}_{2}$ represent the signal and noise parts of $\boldsymbol{x}$. Note that $\boldsymbol{A}_{1}$ and $\boldsymbol{A}_{2}$ are identifiable only up to postmultiplication by $q\times q$ and $(p-q)\times(p-q)$ orthogonal matrices, respectively. If $q\geq p-1$ and the components of $\boldsymbol{z}$ are independent, the model is called independent component model, $\boldsymbol{A}=(\boldsymbol{A}_{1},\boldsymbol{A}_{2})$ is then identified up to the signs and permutation of its columns. Inference on $\boldsymbol{A}$ or its inverse, the unmixing matrix $\boldsymbol{A}^{-1}$, is then known as independent component analysis (ICA). [1] and [2] assumed independent components but allowed $1\leq q\leq p$; we call this approach non-Gaussian independent component analysis (NGICA). In our model there is no restriction on the number of Gaussian components and the non-Gaussian signal components can be dependent of each other. The Gaussian Mixture Models (GMM) with equal covariance matrices and $q+1$ groups and multivariate skew-normal distributions ($q=1$), for example, are included in this wider model. Our definition of the NGCA model is as in [3] and originally suggested in [4]. For recent contributions and references for NGCA but always with known $q$, see e.g. [5] and [6].

In the independent component analysis (ICA) the fourth-order blind identification (FOBI) by [7] uses the regular covariance matrix

[TABLE]

and the scatter matrix based on fourth moments

[TABLE]

where $r^{2}=(\boldsymbol{x}-E(\boldsymbol{x}))^{\prime}\boldsymbol{S}_{1}^{-1}(\boldsymbol{x}-E(\boldsymbol{x}))$, and finds an unmixing matrix $\boldsymbol{W}\in\mathbb{R}^{p\times p}$ such that

[TABLE]

for some diagonal matrix $\boldsymbol{D}=diag(d_{1},...,d_{p})$. $\boldsymbol{W}$ and $\boldsymbol{D}$ then provide the eigenvectors and eigenvalues of $\boldsymbol{S}_{1}^{-1}\boldsymbol{S}_{2}$ and $d_{i}=E[(\boldsymbol{W}\boldsymbol{x})_{i}^{4}]+p-1$, $i=1,...,p$. If in the independent component model the fourth moments of $z_{1},...,z_{p}$ are distinct, $\boldsymbol{W}$ is uniquely defined up to signs and permutations of its rows and $\boldsymbol{W}\boldsymbol{x}$ has independent components with distributions of $z_{1},...,z_{p}$ (up to signs). For a Gaussian $z_{i}$, $d_{i}=p+2$. The matrix $\boldsymbol{D}$ also lists the eigenvalues of the symmetric matrix

[TABLE]

Our test constructions for the dimension of the signal space in the wider NGCA model are based on the estimated eigenvalues of $\boldsymbol{R}$. The test statistics were proposed already in [8] for the NGICA model but without a careful analysis of their limiting distributions.

The plan in this paper is as follows. In Section II our test statistic for testing whether the dimension of the signal space is $k$ is introduced. The asymptotic and bootstrap test versions are provided in Sections III and  IV, respectively. Also the estimate based on the asymptotic test is discussed. The two strategies are compared in simulations in Section V and the paper ends with a discussion on alternative tests. The proofs are provided in the Appendix.

Throughout the paper we use the following notation. The first and second moments of the eigenvalues of positive definite and symmetric $\boldsymbol{R}\in\mathbb{R}^{p\times p}$ are denoted by

[TABLE]

and the variance of the eigenvalues is $s^{2}(\boldsymbol{R}):=m_{2}(\boldsymbol{R})-m_{1}^{2}(\boldsymbol{R})$. We write $\mathcal{O}^{p\times k}$ for the set of $p\times k$ matrices with orthonormal columns, $k\leq p$.

II Test statistic for the dimension

Recall that $\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$ is a random sample from a distribution of $\boldsymbol{x}=\boldsymbol{A}\boldsymbol{z}+\boldsymbol{b}$ where $\boldsymbol{A}\in\mathbb{R}^{p\times p}$ is non-singular, $\boldsymbol{b}\in\mathbb{R}^{p}$, $\mathbb{E}(\boldsymbol{z})=\boldsymbol{0}$ and $Cov(\boldsymbol{z})=\boldsymbol{I}_{p}$. Further $\boldsymbol{z}=(\boldsymbol{z}_{1}^{\prime},\boldsymbol{z}_{2}^{\prime})^{\prime}$ where $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ are independent, $\boldsymbol{z}_{1}\in\mathbb{R}^{q}$ is non-Gaussian and $\boldsymbol{z}_{2}\in\mathbb{R}^{p-q}$ is Gaussian. We also need to assume that the fourth moments exist and $E[(\boldsymbol{u}^{\prime}\boldsymbol{z}_{1})^{4}]\neq 3$ for all $\boldsymbol{u}^{\prime}\boldsymbol{u}=1$, $\boldsymbol{u}\in\mathbb{R}^{q}$. This is a bit stronger assumption than non-Gaussianity of $\boldsymbol{z}_{1}$. With unknown $q$, we then wish to test the null hypothesis

[TABLE]

stating that the dimension of the signal space is $k$.

Natural estimates of $\boldsymbol{S}_{1}$, $\boldsymbol{S}_{2}$ and $\boldsymbol{R}$ are

[TABLE]

[TABLE]

and $\widehat{\boldsymbol{R}}=\widehat{\boldsymbol{S}}_{1}^{-1/2}\widehat{\boldsymbol{S}}_{2}\widehat{\boldsymbol{S}}_{1}^{-1/2}$, respectively. To test the null hypothesis $H_{0,k}$, we use the test statistic

[TABLE]

If

[TABLE]

then $T_{k}=m_{2}\left(\widehat{\boldsymbol{U}}_{k}^{\prime}(\widehat{\boldsymbol{R}}-(p+2)\boldsymbol{I}_{p})\widehat{\boldsymbol{U}}_{k}\right)$.

Note that [9] used $m_{2}\left(\widehat{\boldsymbol{R}}-(p+2)\boldsymbol{I}_{p}\right)$ to test for (full) multivariate normality which is a special case here. Further note that the estimated projections (with respect to Mahalanobis inner product) to the noise and signal subspaces are given by $\widehat{\boldsymbol{Q}}_{k}=\widehat{\boldsymbol{S}}_{1}^{1/2}\widehat{\boldsymbol{U}}_{k}\widehat{\boldsymbol{U}}_{k}^{\prime}\widehat{\boldsymbol{S}}_{1}^{-1/2}$ and $\widehat{\boldsymbol{P}}_{k}=\boldsymbol{I}_{p}-\widehat{\boldsymbol{Q}}_{k}$, respectively.

III Asymptotic test for dimension

As the test statistics $T_{k}$ are invariant under affine transformations $\boldsymbol{A}\boldsymbol{x}_{i}+\boldsymbol{b}$, $i=1,...,n$, we can without loss of generality assume that $\boldsymbol{A}=\boldsymbol{I}_{p}$ and $\boldsymbol{b}=\boldsymbol{0}$. In the following, we consider the limiting behavior of $n(p-k)T_{k}$, $k=0,...,p-1$, under true $H_{0,q}$.

For $k=0,1,...,p-1$, write $\boldsymbol{V}_{k}=(\boldsymbol{0},\boldsymbol{I}_{p-k})^{\prime}\in\mathcal{O}^{p\times(p-k)}$ and

[TABLE]

Then $T_{k}\leq T_{k}^{*}$, $k=0,...,p-1$, and, for the true value $q$, $nT_{q}=nT_{q}^{*}+o_{P}(1)$, see the Appendix. We then have the following.

Theorem 1**.**

*Under the previously stated assumptions and under $H_{0q}$,

(i) for $k<q$, $T_{k}\to_{P}c$ for some $c>0$,

(ii) for $k=q$, $n(p-k)T_{k}\to_{d}C_{k}$, and

(iii) for $k>q$, $n(p-k)T_{k}\leq n(p-k)T_{k}^{*}\to_{d}C_{k}$

where*

[TABLE]

with independent chi squared variables $\chi^{2}_{(p-k-1)(p-k+2)/2}$ and $\chi^{2}_{1}$ and $\sigma_{1}=Var\left(\|\boldsymbol{z}\|^{2}\right)+8$ and $\sigma_{2}=4$.

In the regular testing procedure, the null hypothesis $H_{0,k}$ is the rejected if

[TABLE]

where the critical point $c_{k,\alpha}$ is determined by $\mathbb{P}(C_{k}\geq c_{k,\alpha})=\alpha$. Note that the test, although constructed for $H_{0,k}$ can actually be viewed as a consistent test also for $H^{*}_{0,k}$: At least $p-k$ eigenvalues of $\boldsymbol{R}$ equal $p+2$. To find $c_{k,\alpha}$ in practice, $\sigma_{1}$ must be replaced by its consistent estimate $\hat{\sigma}_{1}$. Then, with increasing $n$, (i) for $k<q$, the power of the test for $H_{0,k}$ goes to one, (ii) for $k=q$, the size of the test for $H_{0,q}$ goes to prespecified $\alpha$, and (iii) for $k>q$ the rejection probability for $H_{0,k}$ tends to be smaller than $\alpha$.

We next discuss the estimation of $\sigma_{1}$. Write $\hat{\boldsymbol{z}}_{i}=\widehat{\boldsymbol{W}}\boldsymbol{x}_{i}$, $i=1,...,n$. Then, even without knowing the true value of $q$, the parameter $\sigma_{1}$ can be consistently estimated by $\hat{\sigma}_{1}=\frac{1}{n}\sum_{i=1}^{n}\|\hat{\boldsymbol{z}}_{i}\|^{4}-p^{2}+8$. Note that in the independent component model, we simply have $\sigma_{1}=\sum_{k=1}^{p}E(z_{k}^{4})-p+8$ which yields as estimate $\hat{\sigma}_{1}=\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{p}(\hat{z}_{i})_{k}^{4}-p+8$. With known $q$, the estimates can be simplified further.

A consistent estimate $\hat{q}$ of the unknown dimension $q$ can be based on the test statistics $T_{k}$, $k=0,...,p-1$ as follows.

Corollary 1**.**

For all $k=0,...,p-1$, let $(c_{k,n})$ be a sequence such that $c_{k,n}\to\infty$ and $\frac{c_{k,n}}{n}\to 0$ as $n\to\infty$. Then

[TABLE]

and

[TABLE]

IV Bootstrap test for dimension

Theorem 1 shows that the limiting distribution of $n(p-q)T_{q}$ (with estimated noise subspace) and $n(p-q)T_{q}^{*}$ (with known noise subspace) are the same. This means that, if one uses the asymptotic test in the small sample case, the variation coming from the estimation of the subspace is ignored. We therefore propose that the small sample null distribution of $n(p-k)T_{k}$ should be estimated by resampling from a distribution for which the null hypothesis $H_{0,k}$ is true and which is as similar as possible to the empirical distribution of $\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$. For this type of bootstrap sampling, we also need estimated projections $\widehat{\boldsymbol{Q}}_{k}=\widehat{\boldsymbol{S}}_{1}^{1/2}\widehat{\boldsymbol{U}}_{k}\widehat{\boldsymbol{U}}_{k}^{\prime}\widehat{\boldsymbol{S}}_{1}^{-1/2}$ and $\widehat{\boldsymbol{P}}_{k}=\boldsymbol{I}_{p}-\widehat{\boldsymbol{Q}}_{k}$. [8] suggested the following procedure for the NGCA model.

Generating a bootstrap sample for $H_{0,k}$:

Starting with centered $\boldsymbol{X}\in\mathbb{R}^{n\times p}$, compute $\widehat{\boldsymbol{S}}_{1}$, $\widehat{\boldsymbol{S}}_{2}$, $\widehat{\boldsymbol{R}}$, $\widehat{\boldsymbol{U}}_{k}$, $\widehat{\boldsymbol{Q}}_{k}$ and $\widehat{\boldsymbol{P}}_{k}$. 2. 2.

Take a bootstrap sample $\widetilde{\boldsymbol{X}}=(\tilde{\boldsymbol{x}}_{1},\ldots,\tilde{\boldsymbol{x}}_{n})^{\prime}$ of size $n$ from $\boldsymbol{X}$. 3. 3.

For the $(p-k)$-dimensional noise space to be gaussian, transform

[TABLE]

and $\boldsymbol{o}_{1},\ldots,\boldsymbol{o}_{n}$ are iid from $N_{p-k}(\boldsymbol{0},\boldsymbol{I}_{p-k})$. 4. 4.

$\boldsymbol{X}^{*}=(\boldsymbol{x}_{1}^{*},...,\boldsymbol{x}_{n}^{*})^{\prime}$.

In [8] bootstrap sampling for the NGICA model was suggested as well.

Let $T=T(\boldsymbol{X})$ be a test statistic for $H_{0,k}$ such as $T_{k}$. If $\boldsymbol{X}_{1}^{*},...,\boldsymbol{X}_{M}^{*}$ are independent bootstrap samples as described above and $T_{i}^{*}=T(\boldsymbol{X}_{i}^{*})$ then the bootstrap $p$-value is given by

[TABLE]

V A simulation study

To compare the asymptotic and bootstrap tests we generate data sets obeying the following three NGCA models. As the test are affine invariant, it is not a restriction to use $\boldsymbol{A}=\boldsymbol{I}_{p}$ and $\boldsymbol{b}=\boldsymbol{0}$ in simulations.

M1:

A GMM model that is a mixture of $N_{6}(6\boldsymbol{e}_{1},\boldsymbol{I}_{6})$, $N_{6}(4\boldsymbol{e}_{2},\boldsymbol{I}_{6})$ and $N_{6}(\boldsymbol{0},\boldsymbol{I}_{6})$ with proportions 0.1, 0.4 and 0.5 respectively. Then $p=6$ and $q=2$.

M2:

An NGCA model with two independent bivariate nongaussian components representing the Greek letters $\boldsymbol{\mu}$ and $\boldsymbol{\Omega}$, see Fig.1, and independent noise $N_{3}(\boldsymbol{0},\boldsymbol{I}_{3})$. Therefore $p=7$ and $q=4$.

M3:

An NGICA model with the non-Gaussian independent components: exponential, $\chi_{1}^{2}$ and uniform and three Gaussian components $N(0,1)$. Hence $p=6$ and $q=3$.

The non-Gaussian components of M1 and M2 are visualized in Figure 1 with samples of size $n=1000$. The complete simulation was performed in R 3.3.2 [10] using the R package ICtest [11] which provides implementations for all methods discussed here. In the comparisons the rejection rates for selected null hypotheses $H_{0,k}$ are reported for asymptotic and bootstrap tests with the test size $\alpha=0.05$. All results are based on 1000 repetitions.

From the simulations we can conclude that for small sample sizes bootstrapping keeps better the target size 0.05 under the true null hypothesis. The sample size needed for decent power naturally depends strongly on the underlying model; M1 seems to be the most difficult case and requires at least 5000 observations. In general, the results are as suggested by the theory.

VI Final remarks

In most applications with high-dimensional data only the non-Gaussian variation is informative and the Gaussian part simply presents noise. Wide literature on NGCA (and NGICA) provides tools to estimate the non-Gaussian subspace but so far always with known dimension $q$. In this paper we suggest efficient asymptotic and bootstrap tests for the dimension that are simply based on the eigenvalues of the easily computable FOBI matrix. Natural estimates of $q$ are also found by successive testing for hypotheses $H_{0,k}$, $k=0,...,p-1$. Note that the stochastic variability of the $p-k$ eigenvectors depends strongly on how close together the corresponding eigenvalues are. In a similar context [12, 13] then suggested a kind of dual estimates of $q$ that are based on the bootstrap variation of eigenvector estimates.

We end the discussion with some remarks on alternative and competing test statistics. The test statistic $T_{k}$ can also be written as a sum $T_{k}=T_{k,1}+T_{k,2}$ where $T_{k,1}=s^{2}(\widehat{\boldsymbol{U}}_{k}^{\prime}\widehat{\boldsymbol{R}}\widehat{\boldsymbol{U}}_{k})$ and $T_{k,2}=[m_{1}(\widehat{\boldsymbol{U}}_{k}^{\prime}\widehat{\boldsymbol{R}}\widehat{\boldsymbol{U}}_{k})-(p+2)]^{2}$. The first part $T_{k,1}$ provides a test statistic for the equality of $p-k$ eigenvalues closest to $p+2$ and the second part $T_{k,2}$ measures the deviation of the average of those eigenvalues from $p+2$ (Gaussian case). Besides $n(p-k)T_{k}$, one can also use

[TABLE]

or their sum as a test statistic. See the Appendix. Under true $H_{0,k}$, these statistics have limiting chi square distributions with $(p-k-1)(p-k+2)/2$, 1, and $(p-k-1)(p-k+2)/2+1$ degrees of freedom. The first two statistics use less information and are therefore in most cases less powerful than their sum, and the behavior of their sum is very similar to that of $T_{k}$ as also seen in our simulations (but not reported here).

Note also that FOBI is just a simple special case of the so called two-scatter method [14, 15]; $\boldsymbol{S}_{1}$ and $\boldsymbol{S}_{2}$ are then replaced by any two scatter matrices specific to the problem at hand. Deriving asymptotic tests or using bootstrap testing strategy with

[TABLE]

for any choices of $\boldsymbol{S}_{1}$ and $\boldsymbol{S}_{2}$ is also possible. The properties of these tests with corresponding estimates is a part of our future work.

VII Appendix: Proof of Theorem 1

For the limiting distributions of the scatter matrices $\widehat{\boldsymbol{S}}_{1}$ and $\widehat{\boldsymbol{S}}_{2}$, we need to assume that the fourth moments of $\boldsymbol{z}_{1}$ exist. Naturally $\boldsymbol{z}_{2}$ has moments of any order. Let $\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$ be a random sample from the distribution of

[TABLE]

where $E(\boldsymbol{z})=\boldsymbol{0}$ and $Cov(\boldsymbol{z})=\boldsymbol{I}_{p}$, $\boldsymbol{z}=(\boldsymbol{z}_{1}^{\prime},\boldsymbol{z}_{2})^{\prime}$, $\boldsymbol{z}_{1}\in\mathbb{R}^{q}$ and $\boldsymbol{z}_{2}\in\mathbb{R}^{p-q}$ are independent, $\boldsymbol{z}_{1}$ is non-Gaussian and $\boldsymbol{z}_{2}$ is Gaussian. Due to affine invariance of the test statistic, it is not a restriction to assume in the following that $\boldsymbol{A}=\boldsymbol{I}_{p}$ and $\boldsymbol{b}=\boldsymbol{0}$.

We write $\widehat{\boldsymbol{S}}_{1}=\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})^{\prime}$ and $\widehat{\boldsymbol{S}}_{2}=\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})^{\prime}\widehat{\boldsymbol{S}}^{-1}(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})(\boldsymbol{x}_{i}-\bar{\boldsymbol{x}})^{\prime}$ and, for known $\boldsymbol{A}=\boldsymbol{I}_{p}$ and $\boldsymbol{b}=\boldsymbol{0}$, we have

[TABLE]

Let

[TABLE]

be partitioned as

[TABLE]

Theorem 1 is then implied by the following four Lemmas, starting with a linearization result for $\widehat{\boldsymbol{R}}_{22}\in\mathbb{R}^{(p-q)\times(p-q)}$.

Lemma 1**.**

Under the stated assumptions,

[TABLE]

where

[TABLE]

Write $\boldsymbol{V}_{k}=(\boldsymbol{0},\boldsymbol{I}_{p-k})^{\prime}\in\mathcal{O}^{p\times(p-k)}$. As, for all $k\geq p$,

[TABLE]

for all $\boldsymbol{U}\in\mathcal{O}^{(p-k)\times(p-k)}$ and $\widetilde{\boldsymbol{R}}$ is the average of iid matrices, we further obtain the following.

Lemma 2**.**

Under the stated assumptions and $k\geq q$,

[TABLE]

where $vec(\boldsymbol{N}_{k})$ has a $(p-k)\times(p-k)$-variate normal distribution with zero mean vector and covariance matrix

[TABLE]

$\boldsymbol{K}_{p,p}=\sum_{i=1}^{p}\sum_{j=1}^{p}(\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{\prime})\otimes(\boldsymbol{e}_{j}\boldsymbol{e}_{i}^{\prime})$* is the commutation matrix, and*

[TABLE]

For $k\geq q$, write

[TABLE]

As in [8], we can show the following

Lemma 3**.**

Under the stated assumptions and $k\geq q$, the random variables $nT_{k,1}^{*}$ and $nT_{k,2}^{*}$ are asymptotically independent and

[TABLE]

and

[TABLE]

where $\sigma_{1}=Var\left(\|\boldsymbol{z}\|^{2}\right)+8$ and $\sigma_{2}=4$.

Lemma 4**.**

Under the stated assumptions, $T_{k}\leq T_{k}^{*}$, for all $k\geq q$, and

[TABLE]

The first part in the Lemma 4 is trivial, the second part follows from Lemma 3.1 in [16].

Acknowledgements

The research of K. Nordhausen, H. Oja and J. Virta was partially supported by the Academy of Finland (grant 268703). D.E. Tyler’s research was partially supported by the National Science Foundation Grant No. DMS-1407751. Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect those of the National Science Foundation.