Discriminant analysis in small and large dimensions

arXiv:1705.02826v1 [] 8 May 2017

Taras Bodnara,1 , Stepan Mazurb , Edward Ngailoa and Nestor Parolyac

a
b

Department of Mathematics, Stockholm University, Roslagsvägen 101, SE-10691 Stockholm, Sweden
Unit of Statistics, School of Business, Örebro University, Fakultetsgatan 1, SE-70182 Örebro, Sweden
c

Institute of Statistics, Leibniz University Hannover, Königsworther Platz 1, D-30167 Hannover,
Germany

Abstract
We study the distributional properties of the linear discriminant function under
the assumption of normality by comparing two groups with the same covariance
matrix but different mean vectors. A stochastic representation for the discriminant function coefficients is derived which is then used to obtain their asymptotic
distribution under the high-dimensional asymptotic regime. We investigate the performance of the classification analysis based on the discriminant function in both
small and large dimensions. A stochastic representation is established which allows
to compute the error rate in an efficient way. We further compare the calculated
error rate with the optimal one obtained under the assumption that the covariance
matrix and the two mean vectors are known. Finally, we present an analytical
expression of the error rate calculated in the high-dimensional asymptotic regime.
The finite-sample properties of the derived theoretical results are assessed via an
extensive Monte Carlo study.

ASM Classification: 62H10, 62E15, 62E20, 60F05, 60B20
Keywords: discriminant function, stochastic representation, large-dimensional asymptotics, random matrix theory, classification analysis

1

Corresponding Author: Taras Bodnar. E-Mail: taras.bodnar@math.su.se. Tel: +46 8 164562. Fax:
+46 8 612 6717. This research was partly supported by the Swedish International Development Cooperation Agency (SIDA) through the UR-Sweden Programme for Research, Higher Education and Institutional Advancement. Stepan Mazur acknowledges financial support from the project ”Ambit fields:
probabilistic properties and statistical inference” funded by Villum Fonden

1

1

Introduction

In the modern world of science and technology, high-dimensional data are present in
various fields such as finance, environment science and social sciences. In the sense of
many complex multivariate dependencies observed in data, formulating correct models and
developing inferential procedures are the major challenges. The traditional multivariate
analysis considers fixed or small sample dimensions, while sample sizes approaching to
infinity. However, its methods cannot longer be used in the high-dimensional setting
where the dimension is not treated as fixed but it is allowed to be comparable to the
sample size.
The covariance matrix is one of the mostly used way to capture the dependence between variables. Although its application is restricted only to linear dependence and
more sophisticated methods, like copula, should be applied in the general case, modeling
dynamics in the covariance matrix is still a very popular subject in both statistics and
econometrics. Recently, a number of papers have been published which deal with estimating the covariance matrix (see, e.g., Ledoit and Wolf (2003), Cai and Liu (2011a), Cai
et al. (2011), Agarwal et al. (2012), Fan et al. (2008), Fan et al. (2013), Bodnar et al.
(2014, 2016)) and testing its structure (see, e.g., Johnstone (2001), Bai et al. (2009), Chen
et al. (2010), Cai and Jiang (2011), Jiang and Yang (2013), Gupta and Bodnar (2014))
in large dimension.
In many applications, the covariance matrix is accompanied by the mean vector. For
example, the product of the inverse sample covariance matrix and the difference of the
sample mean vectors is present in the discriminant function where a linear combination
of variables (discriminant function coefficients) is determined such that the standardized
distance between the groups of observations is maximized. A second example arises in
portfolio theory, where the vector of optimal portfolio weights is proportional to the
products of inverse sample covariance matrix and the sample mean vector (see Bodnar
and Okhrin (2011)).
The discriminant analysis is a multivariate technique concerned with separating distinct sets of objects (or observations) (Johnson et al. (2007)). Its two main tasks are to
distinguish distinct sets of observations and to allocate new observations to previously
defined groups (Rencher and Christensen (2012)). The main methods of the discriminant
analysis are the linear discriminant function and the quadratic discriminant function.
The linear discriminant function is a generalization of Fisher linear discriminant analysis,
a method used in statistics, pattern recognition and machine learning to find a linear
combination of features that characterizes or separates two or more groups of objects
in the best way. The application of the linear discriminant function is restricted to the
assumption of the equal covariance matrix in the groups to be separated. Although the
quadratic discriminant function can be used when the latter assumption is violated, its
application is more computational exhaustive, needs to estimate the covariance matrices
2

of each group, and requires more observations than in the case of linear discriminant function (Narsky and Porter (2013)). Moreover, the decision boundary is easy to understand
and to visualize in high-dimensional settings, if the linear discriminant function is used.
The discriminant analysis is a well established topic in multivariate statistics. Many
asymptotic results are available when the sample sizes of groups to be separated are assumed to be large, while the number of variables is fixed and significantly smaller than
the sample size (see, e.g., Muirhead (1982), Rencher and Christensen (2012)). However, these results cannot automatically be transferred when the number of variables is
comparable to the sample size which is known in the statistical literature as the highdimensional asymptotic regime. It is remarkable that in this case the results obtained
under the standard asymptotic regime can deviate significantly from those obtained under the high-dimensional asymptotics (see, e.g., Bai and Silverstein (2010)). Fujikoshi
and Seo (1997) provided an asymptotic approximation of the linear discriminant function
in high dimension by considering the case of equal sample sizes and compared the results
with the classical asymptotic approximation by Wyman et al. (1990). For the samples of
non-equal sizes, they pointed out that the high-dimensional approximation is extremely
accurate. However, Tamatani (2015) showed that the Fisher linear discriminant function
performs poorly due to diverging spectra in the case of large-dimensional data and small
sample sizes. Bickel and Levina (2004), Srivastava and Kubokawa (2007) investigated
the asymptotic properties of the linear discriminant function in high dimension, while
modifications of the linear discriminant function can be found in Cai and Liu (2011b),
Shao et al. (2011). The asymptotic results for the discriminant function coefficients in
matrix-variate skew models can be found in Bodnar et al. (2017b).
We contribute to the statistical literature by deriving a stochastic representation of
the discriminant function coefficient and the classification rule based on the linear discriminant function. These results provide us an efficient way of simulating these random
quantities and they are also used in the derivation of their high-dimensional asymptotic
distributions, using which the error rate of the classification rule based on the linear discriminant function can be easily assessed and the problem of the increasing dimensionality
can be visualized in a simple way.
The rest of the paper is organized as follows. The finite-sample properties of the discriminant function are presented in Section 2.1, where, in particular we derive a stochastic
representation for the discriminant function coefficients. In Section 2.2, an exact one-sided
test for the comparison of the population discriminant function coefficients is suggested,
while a stochastic representation for the classification rule is obtained in Section 2.3. The
finite-sample results are then use to derive the asymptotic distributions of the discriminant function coefficients and of the classification rule in Section 3, while finite sample
performance of the asymptotic distribution is analysed in Section 3.2.

3

2

Finite-sample properties of the discriminant function
(1)

(1)

(2)

(2)

Let x1 , . . . , xn1 and x1 , . . . , xn2 be two independent samples from the multivariate
normal distributions which consist of independent and identically distributed random
(1)
(2)
vectors with xi ∼ Np (µ1 , Σ) for i = 1, ..., n1 and xj ∼ Np (µ2 , Σ) for j = 1, ..., n2
where Σ is positive definite. Throughout the paper, 1n denotes the n-dimensional vector
of ones, In is the n × n identity matrix, and the symbol ⊗ stands for the Kronecker
product.




(2)
(2)
(1)
(1)
Let X(1) = x1 , . . . , xn1 and X(2) = x1 , . . . , xn2 be observation matrices. Then
the sample estimators for the mean vectors and the covariance matrices constructed from
each sample are given by
x̄

(j)

nj
1
1 X (j)
xi = X(j) 1nj
=
nj i=1
nj
nj

(j)

S

T

1 X  (j)
(j)
.
xi − x̄(j) xi − x̄(j)
=
nj − 1 i=1

The pooled estimator for the covariance matrix, i.e., an estimator for Σ obtained from
two samples, is then given by
Spl =



1
(n1 − 1)S(1) + (n2 − 1)S(2)
n1 + n2 − 2

(1)

The following lemma (see, e.g., (Rencher and Christensen, 2012, Section 5.4.2)) presents
the joint distribution of x̄(1) , x̄(2) and Spl .


Lemma 1. Let X1 ∼ Np,n1 µ1 1Tn1 , Σ ⊗ In1 and X2 ∼ Np,n2 µ2 1Tn2 , Σ ⊗ In2 for p <
n1 + n2 − 2. Assume that X1 and X2 are independent. Then


1
(1)
(a) x̄ ∼ Np µ1 , n1 Σ ,


(b) x̄(2) ∼ Np µ2 , n12 Σ ,
(c) (n1 + n2 − 2)Spl ∼ Wp (n1 + n2 − 2, Σ),
Moreover, x̄(1) , x̄(2) and Spl are mutually independently distributed.
The results of Lemma 1, in particular, implies that
x̄

(1)

− x̄

(2)



 
1
1
∼ Np µ1 − µ2 ,
+
Σ
n1 n2

which is independent of Spl .
4

(2)

2.1

Stochastic representation for the discriminant function coefficients

The discriminant function coefficients are given by the following vector
â = S−1
x̄(1) − x̄(2)
pl



(3)

which is the sample estimator of the population discriminant function coefficient vector
expressed as
a = Σ−1 (µ1 − µ2 )
We consider a more general problem by deriving the distribution of linear combinations
of the discriminant function coefficients. This result possesses several practical application: (i) it allows a direct comparison of the population coefficients in the discriminant
function by deriving a corresponding statistical test; (ii) it can be used in the classification problem where providing a new observation vector one has to decide to which of two
groups the observation vector has to be ordered.
Let L be a k × p matrix of constants such that rank(L) = k < p. We are then
interested in

θ̂ = Lâ = LS−1
x̄(1) − x̄(2) .
pl

(4)

Choosing different matrices L we are able to provide different inferences about the linear
combinations of the discriminant function coefficients. For instance, if k = 1 and L is the
vector with all elements zero except the one on the jth position which is one, then we get
the distribution of the jth coefficient in the discriminant function. If we choose k = 1 and
L = (1, −1, 0, . . . , 0)T , then we analyse the difference between the first two coefficients
in the discriminant function. The corresponding result can be further used to test if the
population counterparts to these coefficients are zero or not. For k > 1 several linear
combinations of the discriminant function coefficients are considered simultaneously.
In the next theorem we derive a stochastic representation for θ̂. The stochastic representation is a very important tool in analysing the distributional properties of random
quantities. It is widely spread in the computation statistics (e.g., Givens and Hoeting
(2012)), in the theory of elliptical distributions (see, Gupta et al. (2013)) as well as in
d
Bayesian statistics (cf., Bodnar et al. (2017a)). Later on, we use the symbol = to denote
the equality in distribution.
Theorem 1. Let L be an arbitrary k × p matrix of constants such that rank(L) = k < p.
Then, under the assumption of Lemma 1 the stochastic representation of θ̂ = Lâ is given

5

by
s
d

θ̂ = (n1 + n2 − 2)ξ −1 LΣ−1 x̆ +

!

x̆T Σ−1 x̆
1/2
LRx̆ LT
t0 ,
n1 + n2 − p



where Rx̆ = Σ−1 −Σ−1 x̆x̆T Σ−1 /x̆T Σ−1 x̆; ξ ∼ χ2n1 +n2 −p−1 , x̆ ∼ Np µ1 − µ2 , n11 +
and t0 ∼ tk (n1 + n2 − p, 0k , Ik ). Moreover, ξ, x̆ and t0 are mutually independent.

(5)

1
n2

 
Σ ,

Proof. From Lemma 1.(c) and Theorem 3.4.1 of Gupta and Nagar (2000) we obtain that
1
S−1 ∼ IW p (n1 + n2 + p − 1, Σ−1 ).
n1 + n2 − 2 pl

(6)

Also, since x̆ = x̄(1) − x̄(2) and Spl are independent, the conditional distribution of
∗
−1 ∗
∗
θ̂ = LS−1
pl x̆ given x̆ = x̆ equals to the distribution of θ = LSpl x̆ and it can be rewritten
in the following form
θ

∗
x̆∗T S−1
pl x̆
.
= (n1 + n2 − 2)x̆ Σ x̆ ∗T −1 ∗
x̆ Spl x̆ (n1 + n2 − 2)x̆∗T Σ−1 x̆∗
d

∗

∗T

−1 ∗

∗
LS−1
pl x̆

Applying Theorem 3.2.12 of Muirhead (1982) we obtain that
x̆∗T Σ−1 x̆∗
ξ = (n1 + n2 − 2) ∗T −1 ∗ ∼ χ2n1 +n2 −p−1
x̆ Spl x̆
∗

(7)

and its distribution is independent of x̆∗ . Hence,
ξ = (n1 + n2 − 2)

x̆T Σ−1 x̆
∼ χ2n1 +n2 −p−1
x̆T S−1
x̆
pl

(8)

and ξ, x̆ are independent.
∗
Using Theorem 3 of Bodnar and Okhrin (2008) we get that x̆∗T S−1
pl x̆ is indepen∗
∗T −1 ∗
∗
∗
∗T −1 ∗
dent of LS−1
pl x̆ /x̆ Spl x̆ for given x̆ . Therefore, ξ is independent of x̆ Σ x̆ ·
−1
T −1
T −1
∗
∗T −1 ∗
LS−1
pl x̆ /x̆ Spl x̆ and, respectively, ξ is independent of x̆ Σ x̆ · LSpl x̆/x̆ Spl x̆. Furthermore, from the proof of Theorem 1 of Bodnar and Schmid (2008) it holds that
∗
LS−1
pl x̆



∗T −1 ∗
−1 ∗ x̆ Σ x̆
T
x̆ Σ x̆ ∗T −1 ∗ ∼ tk n1 + n2 − p; LΣ x̆ ,
LRx̆∗ L
n1 + n2 − p
x̆ Spl x̆
∗T

−1 ∗

(9)

with Rx̆∗ = Σ−1 − Σ−1 x̆∗ x̆∗T Σ−1 /x̆∗T Σ−1 x̆∗ .
Thus, we obtain the following stochastic representation of θ̂ which is given by
s
d

θ̂ = (n1 + n2 − 2)ξ −1 LΣ−1 x̆ +

6

!

x̆T Σ−1 x̆
1/2
LRx̆ LT
t0 ,
n1 + n2 − p

(10)

 


1
1
where Rx̆ = Σ −Σ x̆x̆ Σ /x̆ Σ x̆; ξ ∼
x̆ ∼ Np µ1 − µ2 , n1 + n2 Σ ,
and t0 ∼ tk (n1 + n2 − p, 0k , Ik ). Moreover, ξ, x̆ and t0 are mutually independent. The
theorem is proved.
−1

−1

T

−1

−1

T

χ2n1 +n2 −p−1 ,

In the next corollary we consider the special case when k = 1, that is, when L = lT is
a p-dimensional vector of constants.
Corollary 1. Let λ = 1/n1 +1/n2 and let l be a p-dimensional vector of constants. Then,
under the condition of Theorem 1, the stochastic representation of θ̂ = lT â is given by
d

θ̂ = (n1 + n2 − 2)ξ −1

s
!

λ(p − 1)
λ+
u lT Σ−1 lz0 , (11)
lT Σ−1 (µ1 − µ2 ) +
n1 + n2 − p


where ξ ∼ χ2n1 +n2 −p−1 , z0 ∼ N (0, 1), u ∼ F p − 1, n1 + n2 − p, (µ1 − µ2 )T Rl (µ1 − µ2 )/λ
(non-central F-distribution with p−1 and n1 +n2 −p degrees of freedom and non-centrality
parameter (µ1 − µ2 )T Rl (µ1 − µ2 )/λ) with Rl = Σ−1 − Σ−1 llT Σ−1 /lT Σ−1 l; ξ, z0 and u
are mutually independently distributed.
Proof. From Theorem 1 we get that
s
d

θ̂ = (n1 + n2 − 2)ξ −1

lT Σ−1 x̆ + t0

x̆T Σ−1 x̆

!

· lT Rx̆ l
n1 + n2 − p


p
√
t0
−1
T −1
= (n1 + n2 − 2)ξ
lT Σ−1 l x̆T Rl x̆ ,
l Σ x̆ + √
n1 + n2 − p

(12)
(13)

where Rl = Σ−1 − Σ−1 llT Σ−1 /lT Σ−1 l; ξ ∼ χ2n1 +n2 −p−1 , t0 ∼ t(n1 + n2 − p, 0, 1), and
x̆ ∼ Np (µ1 − µ2 , λΣ) with λ = 1/n1 + 1/n2 ; ξ, t0 and x̆ are mutually independent.
Because x̆ ∼ Np (µ1 − µ2 , λΣ), Rl ΣRl = Rl , and tr [Rl Σ] = p − 1, the application of
Corollary 5.1.3a of Mathai and Provost (1992) leads to
ζ = λ−1 x̆T Rl x̆ ∼ χ2p−1 δ 2



(14)

where δ 2 = (µ1 − µ2 )T Rl (µ1 − µ2 )/λ. Moreover, since Rl ΣΣ−1 l = 0, the application of
Theorem 5.5.1 of Mathai and Provost (1992) proves that lT Σ−1 x̆ and ζ are independently
distributed.
Finally, we note that the random variable t0 ∼ t(n1 + n2 − p, 0, 1) has the following
stochastic representation
d

t0 =z0

r

n1 + n2 − p
,
w

7

(15)

where z0 ∼ N (0, 1) and w ∼ χ2n1 +n2 −p ; z0 and w are independent. Hence,
s
−1

T

l Σ x̆ + t0

λζ · lT Σ−1 l
ζ, w ∼ N
n1 + n2 − p





ζ
l Σ µ, λl Σ l 1 +
(16)
w



p−1
T −1
T −1
u
,(17)
= N l Σ µ, λl Σ l 1 +
n1 + n2 − p
T

−1

T

−1

where
u=


ζ/(p − 1)
∼ F p − 1, n1 + n2 − p, (µ1 − µ2 )T Rl (µ1 − µ2 )/λ .
w/(n1 + n2 − p)

(18)

Putting all above together we get the statement of the corollary.

2.2

Test for the population discriminant function coefficients

One of the most important questions when the discriminant analysis is performed is
to decide which coefficients are the most influential in the decision. Several methods
exist in the literature with the following three approaches to be the most popular (c.f.,
(Rencher and Christensen, 2012, Section 5.5)): (i) standardized coefficients; (ii) partial
F -values; (iii) correlations between the variables and the discriminant function. (Rencher,
1998, Theorem 5.7A) argued that each of this three methods has several drawbacks. For
instance, the correlations between the variables and the discriminant function do not show
the multivariate contribution of each variable, but provide only univariate information how
each variable separates the groups, ignoring the presence of other variables.
In this section, we propose an alternative approach based on the statistical hypothesis
test. Namely, exact statistical tests will be derived on the null hypothesis that two
population discriminant function coefficients are equal (two-sided test) as well as on the
alternative hypothesis that a coefficient in the discriminant function is larger than another
one (one-sided test). The testing hypothesis for the equality of the i-th and the j-th
coefficients in the population discriminant function is given by
H0 : ai = aj

against H1 : ai 6= aj ,

(19)

while in the case of one-sided test we check if
H0 : ai ≤ aj

against H1 : ai > aj .

8

(20)

In both cases the following test statistic is suggested
(1)
p
lT S−1
− x̄(2) )
pl (x̄
q
T = n1 + n2 − p − 1 q
lT S−1
l
(n1 + n2 − 2)( n11 + n12 ) + (x̄(1) − x̄(2) )T R̂l (x̄(1) − x̄(2) )
pl
(21)
with
> −1
S−1
pl ll Spl
−1
R̂l = Spl − > −1
and l = (0, .., 0, |{z}
1 , 0..., 0, |{z}
−1 , 0, ..., 0)> .
l S l
pl

i

j

The distribution of T follows from (Bodnar and Okhrin, 2011, Theorem 6) and it is
summarized in Theorem 2.
Theorem 2. Let λ = 1/n1 + 1/n2 and let l be a p-dimensional vector of constants. Then,
under the condition of Theorem 1,
(a) the density of T is given by
n1 + n2 − p
fT (x) =
λ(p − 1)

Z
0

∞


ftn1 +n2 −p−1,δ1 (y) (x)fFp−1,n1 +n2 −p,s/λ


n1 + n2 − p
y dy (22)
λ(p − 1)

√
lT Σ−1 (µ1 −µ2 )
√
with δ1 (y) = η/ λ + y, η =
, and s = (µ1 −µ2 )T Rl (µ1 −µ2 ); the symbol
lT Σ−1 l
fG (.) denotes the density of the distribution G.
(b) Under the null hypothesis it holds that T ∼ tn1 +n2 −p−1 and T is independent of (x̄(1) −
x̄(2) )T R̂l (x̄(1) − x̄(2) ).
Theorem 2 shows that the test statistics T has a standard t-distribution under the null
hypothesis. As a result, the suggested test will reject the null hypothesis of the two-sided
test (19) as soon as |T | > tn1 +n2 −p−1;1−α/2 .
The situation is more complicated in the case of the one-sided test (20). In this case
the maximal probability of the type I error has to be control. For that reason, we first
calculate the probability of rejection of the null hypothesis for all possible parameter
values and after that we calculate its maximum for the parameters which correspond to
the null hypothesis in (20). Since the distribution of T depends on µ1 , µ2 , and Σ only
over η and s (see, Theorem 2), the task of finding the maximum is significantly simplified.
Let FG (.) denotes the distribution function of the distribution G. For any constant q, we

9

get
Z

+∞

fT (x)dx

P(T > q) =
q

Z
n1 + n2 − p
n1 + n2 − p ∞
ftn1 +n2 −p−1,δ1 (y) (x)fFp−1,n1 +n2 −p,s/λ (
y)dydx
λ(p − 1) 0
λ(p − 1)
q
 Z +∞

Z
n1 + n2 − p ∞
n1 + n2 − p
y
fFp−1,n1 +n2 −p,s/λ
ftn1 +n2 −p−1,δ1 (y) (x)dxdy
λ(p − 1) 0
λ(p − 1)
q


Z
n1 + n2 − p
n1 + n2 − p ∞
(1 − Ftn1 +n2 −p−1,δ1 (y) (q))fFp−1,n1 +n2 −p,s/λ
y dy
λ(p − 1) 0
λ(p − 1)


Z
n1 + n2 − p ∞
n1 + n2 − p
(1 − Ftn1 +n2 −p−1,0 (q))fFp−1,n1 +n2 −p,s/λ
y dy
λ(p − 1) 0
λ(p − 1)
(1 − Ftn1 +n2 −p−1,0 (q)).
Z

=
=
=
≤
=

+∞

where the last equality follows from the fact that the distribution function of the noncentral t-distribution is a decreasing function in non-centrality parameter and δ1 (y) ≤ 0.
Consequently, we get q = tn1 +n2 −p−1;1−α and the one-sided test rejects the null hypothesis
in (20) as soon as T > tn1 +n2 −p−1;1−α .

2.3

Classification analysis

Having a new observation vector x, we classify it to one of the considered two groups.
Assuming that no prior information is available about the classification result, i.e. the
prior probability of each group is 1/2, the decision which is based on the optimal rule is
to assign the observation vector x to the first group as soon as the following inequality
holds (c.f., (Rencher, 1998, Section 6.2))
1
(µ1 − µ2 )> Σ−1 x > (µ1 − µ2 )> Σ−1 (µ1 + µ2 )
2

(23)

and to the second group otherwise. The error rate is defined as the probability of classifying the observation x into one group, while it comes from another one. Rencher (1998)
presented the expression of the error rate expressed as
1
P(classify to the first group | second group is true)
2
1
+
P(classify to the second group | first group is true)
2

∆
= Φ −
with ∆2 = (µ1 − µ2 )> Σ−1 (µ1 − µ2 ) ,
2

ERp (∆) =

where Φ(.) denotes the distribution function of the standard normal distribution.
In practice, however, µ1 , µ2 , and Σ are unknown quantities and the decision is based

10

on the inequality
1 (1)
(1)
(x̄(1) − x̄(2) )> S−1
− x̄(2) )> S−1
+ x̄(2) )
pl x > (x̄
pl (x̄
2

(24)

instead. Next, we derive the error rate of the decision rule (24). Let
1 (1)
(1)
− x̄(2) )> S−1
+ x̄(2) )
dˆ = (x̄(1) − x̄(2) )> S−1
pl (x̄
pl x − (x̄
2


1 (1)
(1)
(2) > −1
(2)
= (x̄ − x̄ ) Spl x − (x̄ + x̄ ) .
2

(25)

ˆ
In Theorem 3 we present the stochastic representation of d.
Theorem 3. Let λ = 1/n1 +1/n2 . Then, under the condition of Theorem 1, the stochastic
representation of dˆ is given by

 (−1)i−1 
√
√
n1 + n2 − 2
λni − 2 
d
dˆ =
(−1)i−1
λξ2 + (∆ + λw0 )2 +
∆2 + λ∆w0
ξ
2λni
λni
s
!
q
√
p−1
1
+
+
u
1+
λξ2 + (∆ + λw0 )2 z0
for i = 1, 2, (26)
n1 + n2 n1 + n2 − p
where u|ξ1 , ξ2 , w0 ∼ F (p − 1, n1 + n2 − p, (n1 + n2 )−1 ξ1 ) with ξ1 |ξ2 , w0 ∼ χp−1,δξ2

2 ,w0

∆2 ξ√
2

and

δξ22 ,w0 = nn1 n2 2 λξ +(∆+ λw )2 , z0 , w0 ∼ N (0, 1), ξ ∼ χ2n1 +n2 −p−1 , ξ2 ∼ χ2p−1 ; ξ, z0 are inde2
0
i
pendent of u, ξ1 , ξ2 , w0 where ξ2 and w0 are independent as well.
Proof. Let x ∼ Np (µi , Σ), Since x̄(1) , x̄(2) , x, and Spl are independently distributed, we
(1)
(2)
get that the conditional distribution of dˆ given x̄(1) = x0 and x̄(2) = x0 is equal to the
distribution of d0 defined by
(1)

(2)

d0 = (x̄0 − x̄0 )> S−1
pl x̃ ,


(1)
(2)
(1)
(2)
where x̃ = x − 12 (x̄0 + x̄0 ) ∼ Np µi − 21 (x̄0 + x̄0 ), Σ , (n1 + n2 − 2)Spl ∼ Wp (n1 +
n2 − 2, Σ), x̃ and Spl are independent.
Following the proof of Corollary 1, we get
d

d0 = (n1 + n2 − 2)ξ
s
+
1+

−1

(1)
(x̄0

−

(2)
x̄0 )T Σ−1



1 (1)
(2)
µi − (x̄0 + x̄0 )
2
!


(p − 1)
(1)
(2)
(1)
(2)
u (x̄0 − x̄0 )T Σ−1 (x̄0 − x̄0 )z0 ,
n1 + n2 − p



T


(1)
(2)
(1)
(2)
1
1
where u ∼ F p − 1, n1 + n2 − p, µi − 2 (x̄0 + x̄0 ) R0 µi − 2 (x̄0 + x̄0 )
with
11

(1)

(2)

(1)

(2)

(1)

(2)

(1)

(2)

R0 = Σ−1 − Σ−1 (x̄0 − x̄0 )(x̄0 − x̄0 )T Σ−1 /(x̄0 − x̄0 )T Σ−1 (x̄0 − x̄0 ), z0 ∼ N (0, 1),
and ξ ∼ χ2n1 +n2 −p−1 ; ξ, z0 and u are mutually independently distributed.
In using that
1 (1)
1 (1)
(2)
(i)
(2)
µi − (x̄0 + x̄0 ) = µi − x̄0 + (−1)i−1 (x̄0 − x̄0 )
2
2
(1)

(2)

and (x̄0 − x̄0 )T R0 = 0, we get

n1 + n2 − 2 (−1)i−1 (1)
d
(x̄ − x̄(2) )T Σ−1 (x̄(1) − x̄(2) ) − (x̄(1) − x̄(2) )T Σ−1 x̄(i) − µi
dˆ =
ξ
2
s
!

p−1
+
1+
u (x̄(1) − x̄(2) )T Σ−1 (x̄(1) − x̄(2) )z0 ,
n1 + n2 − p

T

where u|x̄(1) , x̄(2) ∼ F p − 1, n1 + n2 − p, x̄(i) − µi Rx x̄(i) − µi with Rx = Σ−1 −
Σ−1 (x̄(1) − x̄(2) )(x̄(1) − x̄(2) )T Σ−1 /(x̄(1) − x̄(2) )T Σ−1 (x̄(1) − x̄(2) ), z0 ∼ N (0, 1), and ξ ∼
χ2n1 +n2 −p−1 ; ξ, z0 are independent of u, x̄(1) , x̄(2) .
Since x̄(1) and x̄(2) are independent and normally distributed, we get that
(i)

x̄ − µi
x̄(1) − x̄(2)

!

0
µ1 − µ2

∼ N2p

!
,

1
Σ
ni
i−1
(−1)
Σ
ni

(−1)i−1
Σ
ni

!!

λΣ

and, consequently,
x̄

(i)

− µi |(x̄

(1)

(2)

− x̄ ) ∼ Np




1
(−1)i−1 (1)
(2)
(x̄ − x̄ − (µ1 − µ2 )),
Σ ,
λni
n1 + n2

1
where we used that n1i − λn1 2 = n1 +n
.
2
i
The application of Theorem 5.5.1 in Mathai and Provost (1992) shows that given
(1)
(x̄ − x̄(2) ) the random variables (x̄(1) − x̄(2) )T Σ−1 (x̄(i) − µi ) and (x̄(i) − µi )Rx (x̄(i) − µi )
are independently distributed with

(x̄(1) − x̄(2) )T Σ−1 (x̄(i) − µi )|(x̄(1) − x̄(2) )


(−1)i−1 (1)
1
(2) T −1
(1)
(2)
(1)
(2) T −1
(1)
(2)
∼ N
(x̄ − x̄ ) Σ (x̄ − x̄ − (µ1 − µ2 )),
(x̄ − x̄ ) Σ (x̄ − x̄ )
λni
n1 + n2
and, by using Corollary 5.1.3a of Mathai and Provost (1992),
(n1 + n2 )(x̄(i) − µi )T Rx (x̄(i) − µi )|(x̄(1) − x̄(2) ) ∼ χp−1,δx2

12

with
n1 + n2 (1)
(x̄ − x̄(2) − (µ1 − µ2 ))T Rx (x̄(1) − x̄(2) − (µ1 − µ2 ))
λ2 n2i
n1 + n2
=
(µ1 − µ2 )T Rx (µ1 − µ2 )
λ2 n2i
n1 + n2 (µ1 − µ2 )T Σ−1 (µ1 − µ2 )
(x̄(1) − x̄(2) )T Rµ (x̄(1) − x̄(2) )
=
λ2 n2i (x̄(1) − x̄(2) )T Σ−1 (x̄(1) − x̄(2) )

δx2 =

where we use that (x̄(1) −x̄(2) )T Rx = 0 and Rµ = Σ−1 −Σ−1 (µ1 −µ2 )(µ1 −µ2 )T Σ−1 /(µ1 −
µ2 )T Σ−1 (µ1 − µ2 ).
As a result, we get
n1 + n2 − 2
λni − 2 2 (−1)i−1
d
(−1)i−1
∆x +
(µ1 − µ2 )T Σ−1 (x̄(1) − x̄(2) )
dˆ =
ξ
2λni
λni
s
!

1
p−1
1+
u ∆x z0 ,
+
+
n1 + n2 n1 + n2 − p
where ∆2x = (x̄(1) −x̄(2) )T Σ−1 (x̄(1) −x̄(2) ), u|x̄(1) , x̄(2) ∼ F (p − 1, n1 + n2 − p, (n1 + n2 )−1 ξ1 )
with ξ1 ∼ χp−1,δx2 , z0 ∼ N (0, 1), and ξ ∼ χ2n1 +n2 −p−1 ; ξ, z0 are independent of u, ξ1 , x̄(1) , x̄(2) .
Finally, it holds with ∆2 = (µ1 − µ2 )T Σ−1 (µ1 − µ2 ) that
∆2x

= (x̄

(1)

(2) T

− x̄ ) Rµ (x̄

(1)

2
(µ1 − µ2 )T Σ−1 (x̄(1) − x̄(2) )
− x̄ ) +
,
∆2
(2)

where both summands are independent following Theorem 5.5.1 in Mathai and Provost
(1992). The application of Corollary 5.1.3a in Mathai and Provost (1992) leads to
λ−1 (x̄(1) − x̄(2) )T Rµ (x̄(1) − x̄(2) ) ∼ χ2p−1
and
(µ1 − µ2 )T Σ−1 (x̄(1) − x̄(2) ) ∼ N (∆2 , λ∆2 ).
From the last statement we get the stochastic representation of dˆ expressed as
 (−1)i−1 

√
√
n1 + n2 − 2
λni − 2 
d
dˆ =
(−1)i−1
λξ2 + (∆ + λw0 )2 +
∆2 + λ∆w0
ξ
2λni
λni
s
!
q
√
1
p−1
+
1+
+
u
λξ2 + (∆ + λw0 )2 z0 ,
n1 + n2 n1 + n2 − p
where u|ξ1 , ξ2 , w0 ∼ F (p − 1, n1 + n2 − p, (n1 + n2 )−1 ξ1 ) with ξ1 |ξ2 , w0 ∼ χp−1,δξ2
δξ22 ,w0 =

∆2

n1 +n2
√
λξ2 ,
λ2 n2i λξ2 +(∆+ λw0 )2

2 ,w0

and

z0 , w0 ∼ N (0, 1), ξ ∼ χ2n1 +n2 −p−1 , ξ2 ∼ χ2p−1 ; ξ, z0 are

13

independent of u, ξ1 , ξ2 , w0 where ξ2 and w0 are independent as well.
Theorem 3 shows that the distribution of dˆ is determined by six random variables
ξ, ξ1 , ξ2 , z0 , w0 , and u. Moreover, it depends on µ1 , µ2 , and Σ only via the quadratic form
∆. As a result, the the error rate based on the decision rule (24) is a function of ∆ only
and it is calculated by
ERs (∆) =

1
1 ˆ
P(d > 0| second group is true) + P(dˆ ≤ 0| first group is true) .
2
2

The two probabilities in (27) can easily be approximated for all ∆, p, n1 , and n2 with
high precision by applying the results of Theorem 3 via the following simulation study
(i) Fix ∆ and i ∈ {1, 2}.
(ii) Generate four independent random variables ξb ∼ χ2n1 +n2 −p−1 , ξ2;b ∼ χ2p−1 , z0;b ∼
N (0, 1), and w0;b ∼ N (0, 1).
(iii) Generate ξ1,b ∼ χp−1,δξ2

2 ,w0

with δξ22,b ,w0,b =

∆2 ξ2;b
n1 n2
√
.
n2i λξ2;b +(∆+ λw0;b )2

(iv) Generate u ∼ F (p − 1, n1 + n2 − p, (n1 + n2 )−1 ξ1,b ).
(i)
(v) Calculate dˆb following the stochastic representation (26) of Theorem 3.
(i)
(i)
(vi) Repeat steps (ii)-(v) for b = 1, ..., B leading to the sample dˆ1 , ..., dˆB .

The procedure has to be performed for both values of i = 1, 2 where for i = 1 the relative
number of events {dˆ > 0} will approximate the first summand in (27) while for i = 2 the
relative number of events {dˆ ≤ 0} will approximate the second summand in (27).
It is important to note that the difference between the error rates calculated for the
two decision rules (23) ad (24) could be very large as shown in Figure 1 where ERp (∆) and
ERs (∆) calculated for several values of n1 = n2 ∈ {50, 100, 150, 250} with fixed values
of p ∈ {10, 25, 50, 75}. If p = 10 we do not observe large differences between ERp (∆)
and ERs (∆) computed for different sample sizes. However, this statement does not hold
any longer when p becomes comparable to both n1 and n2 as documented for p = 50 and
p = 75. This case is known in the literature as a large-dimensional asymptotic regime
and it is investigated in detail in Section 3.

3

Discriminant analysis under large-dimensional asymptotics

In this section we derive the asymptotic distribution of the discriminant function coefficients under the high-dimensional asymptotic regime, that is, when the dimension

14

0.5

p=25

0.5

p=10

0.1

0.2

0.3

0.4

ERp(Δ)
ERs(Δ), n1 = n2 = 50
ERs(Δ), n1 = n2 = 100
ERs(Δ), n1 = n2 = 150
ERs(Δ), n1 = n2 = 250

0.0

0.0

0.1

0.2

0.3

0.4

ERp(Δ)
ERs(Δ), n1 = n2 = 50
ERs(Δ), n1 = n2 = 100
ERs(Δ), n1 = n2 = 150
ERs(Δ), n1 = n2 = 250

0

1

2

3

4

5

0

1

2

3

4

5

p=75

0.5

0.5

p=50

0.4

ERp(Δ)
ERs(Δ), n1 = n2 = 50
ERs(Δ), n1 = n2 = 100
ERs(Δ), n1 = n2 = 150
ERs(Δ), n1 = n2 = 250

0.0

0.0

0.1

0.1

0.2

0.2

0.3

0.3

0.4

ERp(Δ)
ERs(Δ), n1 = n2 = 50
ERs(Δ), n1 = n2 = 100
ERs(Δ), n1 = n2 = 150
ERs(Δ), n1 = n2 = 250

0

1

2

3

4

5

0

1

2

3

4

Figure 1: Error rates ERp (∆) and ERs (∆) as functions of ∆ for p ∈ {10, 25, 50, 75} and
ERs (∆).

15

5

increases together with the sample sizes and they all tend to infinity. More precisely, we
assume that p/(n1 + n2 ) → c ∈ [0, 1) as n1 + n2 → ∞.
The following conditions are needed for the validity of the asymptotic results:
(A1) There exists γ ≥ 0 such that p−γ (µ1 − µ2 )T Σ−1 (µ1 − µ2 ) < ∞ uniformly on p.
(A2) 0 <

lim
(n1 ,n2 )→∞

(n1 /n2 ) < ∞.

It is remarkable that, no assumption on the eigenvalues of the covariance matrix Σ,
like they are uniformly bounded on p, is imposed. The asymptotic results are also valid
when Σ possesses unbounded spectrum as well as when its smallest eigenvalue tends to
zero as p → ∞. The constant γ is a technical one and it controls the growth rate of the
quadratic form. In Theorem 4 the asymptotic distribution of linear combinations of the
discriminant function coefficients is provided.
Theorem 4. Assume (A1) and (A2). Let l be a p-dimensional vector of constants such
that p−γ lT Σ−1 l < ∞ is uniformly on p, γ ≥ 0. Then, under the conditions of Theorem 1,
the asymptotic distribution of θ̂ = lT â is given by
√

n1 +

n2 σγ−1


θ̂ −


1 T −1
D
l Σ (µ1 − µ2 ) −→ N (0, 1)
1−c

for p/(n1 + n2 ) → c ∈ [0, 1) as n1 + n2 → ∞ with
σγ2 =

1
(1 − c)3

lT Σ−1 (µ1 − µ2 )

2

+ lT Σ−1 l(µ1 − µ2 )T Σ−1 (µ1 − µ2 )

(27)

!
+ λ(n1 + n2 )lT Σ−1 l1{0} (γ)
where 1A (.) denotes the indicator function of set A.
Proof. Using the stochastic representation (11) of Corollary 1, we get
√


1 T −1
n1 +
l Σ (µ1 − µ2 )
1−c

 −γ T −1
√
1
p l Σ (µ1 − µ2 )
d
−1
=
n1 + n2 (n1 + n2 − 2)ξ −
1−c
p−γ σγ
s
 p −γ T −1
p
n1 + n2 − 2
p
−
1
p l Σ l
+
λ(n1 + n2 )
p−γ + p−γ
u
z0 ,
ξ
n1 + n2 − p
p−γ σγ
n2 σγ−1


θ̂ −


where ξ ∼ χ2n1 +n2 −p−1 , z0 ∼ N (0, 1), u ∼ F p − 1, n1 + n2 − p, (µ1 − µ2 )T Rl (µ1 − µ2 )/λ
with Rl = Σ−1 − Σ−1 llT Σ−1 /lT Σ−1 l; ξ, z0 and u are mutually independently distributed.

16

Since, ξ ∼ χ2n1 +n2 −p−1 , we get that
p
n1 + n2 − p − 1




ξ
D
− 1 −→ N (0, 2)
n1 + n2 − p − 1

for p/(n1 + n2 ) → c ∈ [0, 1) as n1 + n2 → ∞ and, consequently,


√
1
n1 + n2 − p − 1 1
n1 + n2
−1
n1 + n2 (n1 + n2 − 2)ξ −
=√
1−c
ξ
1−c
n1 + n2 − p − 1




p
n1 + n2 − 2
ξ
2
D
×
n1 + n2 − p − 1 (1 − c)
−
−→ z˜0 ∼ N 0,
n1 + n2 − p − 1 n1 + n2 − p − 1
1−c
√

for

p

−γ

p
n1 +n2

= c + o((n1 + n2 )−1/2 ) where z0 and z̃0 are independent.
Furthermore, we get (see, (Bodnar and Reiß, 2016, Lemma 3))
+p

−γ

p−1
c
u − 1{0} (γ) −
n1 + n2 − p
1−c



p−γ (µ1 − µ2 )T Rl (µ1 − µ2 )
1{0} (γ) +
cλ(n1 + n2 )



a.s.

−→ 0

Putting the above results together, we get the statement of the theorem with
σγ2 =

2
1
T −1
2
l
Σ
(µ
−
µ
)
+ lT Σ−1 l(µ1 − µ2 )T Rl (µ1 − µ2 )
1
2
(1 − c)3
!

+ λ(n1 + n2 )lT Σ−1 l1{0} (γ)
=

1
(1 − c)3

lT Σ−1 (µ1 − µ2 )

2

+ lT Σ−1 l(µ1 − µ2 )T Σ−1 (µ1 − µ2 )

!
T

−1

+ λ(n1 + n2 )l Σ l1{0} (γ)

The results of Theorem 4 show that the quantity γ is present only in the asymptotic
variance σγ2 . Moreover, if γ > 0, then the factor λ(n1 + n2 ) vanishes and therefore the
assumption (A2) is no longer needed. However, in the case of γ = 0 we need (A2) in order
to keep the variance bounded. We further investigate this point via simulations in Section
3.3, by choosing γ > 0 and considering small n1 and large n2 such that n1 /n2 → 0.

17

3.1

Classification analysis in high dimension

The error rate of the classification analysis based on the optimal decision rule (23) remains
the same independently of p and it is always equal to


∆
ERp (∆) = Φ −
2



with ∆2 = (µ1 − µ2 )> Σ−1 (µ1 − µ2 ) .

In practice, however, µ1 , µ2 , and Σ are not known and, consequently, one has to make
the decision based on (24) instead of (23). In Theorem 5, we derived the asymptotic
distribution of dˆ under the large-dimensional asymptotics.
e 2 and λni → bi for p/(n1 + n2 ) →
Theorem 5. Assume (A1) and (A2). Let p−γ ∆2 → ∆
c ∈ [0, 1) as n1 + n2 → ∞. Then, under the conditions of Theorem 1, it holds that
pmin(γ,1)/2
D

dˆ
n1 + n2 − 2 (−1)i−1 −γ 2
−
p ∆
pγ
n1 + n2 − p − 1
2

−→ N (−1)i−1

!

c bi − 2
(b1 + b2 )1{0} (γ),
1 − c 2bi

!
c
1
e 4 1[1,+∞) (γ) +
e 2 1[0,1] (γ))
∆
(c(b1 + b2 )1{0} (γ) + ∆
2(1 − c)3
(1 − c)3
for p/(n1 + n2 ) → c ∈ [0, 1) as n1 + n2 → ∞.
Proof. The application of Theorem 3 leads to
!
i−1
ˆ
n
+
n
−
2
(−1)
d
1
2
−
p−γ ∆2
pmin(γ,1)/2
pγ
n1 + n2 − p − 1
2


√
p
n1 + n2 − 2 p
ξ
d
min(γ,1)/2−1/2
√
= p
n1 + n2 − p − 1 1 −
ξ
n1 + n2 − p − 1
n1 + n2 − p − 1
×
×
+

+
×

(−1)i−1 −γ 2 n1 + n2 − 2
λni − 2
p ∆ +
(−1)i−1
2
ξ
2λni


p
√
min(γ,1)/2−γ
min(γ,1)/2−γ/2
min(γ,1)/2−γ
2
−γ
2
p ∆ λw0 + p
λw0
p
λξ2 + 2p
!
(−1)i−1 min(γ,1)/2−γ/2 p −γ 2 √
p
p ∆ λw0
λni
s

n1 + n2 − 2
1
p−1
1+
+
u
ξ
n1 + n2 n1 + n2 − p
!
q
p
√
pmin(γ,1)−2γ λξ2 + (pmin(γ,1)/2−γ/2 p−γ ∆2 + pmin(γ,1)/2−γ λw0 )2 z0

18

D

−→ N (−1)i−1

c bi − 2
(b1 + b2 )1{0} (γ),
1 − c 2bi

!
c
1
e 4 1[1,+∞) (γ) +
e 2 1[0,1] (γ)) ,
∆
(c(b1 + b2 )1{0} (γ) + ∆
2(1 − c)3
(1 − c)3
where the last line follows from Lemma 3 in Bodnar and Reiß (2016) and Slutsky Theorem
(see, (DasGupta, 2008, Theorem 1.5)).
The parameters of the limit distribution derived in Theorem 5 can be significantly
simplified in the special case of n1 = n2 because of λn1 = λn2 = 2. The results of
Theorem 5 are also used to derived the approximate error rate for the decision (24). Let
1 1 −γ
p ∆. Then, the error rate is given by
a = 1−c
2
o 1 n
o
1 nˆ
P d > 0|i = 2 + P dˆ ≤ 0|i = 1
2 (
2
!
)
ˆ
1
d
=
P pmin(γ,1)/2
− (−1)i−1 a > −pmin(γ,1)/2 (−1)i−1 a|i = 2
2
pγ
(
!
)
ˆ
1
d
+
P pmin(γ,1)/2
− (−1)i−1 a ≤ −pmin(γ,1)/2 (−1)i−1 a|i = 1
2
pγ

 min(γ,1)/2



1
ap
− m2
1
−apmin(γ,1)/2 − m1
1−Φ
+ Φ
,
≈
2
v
2
v

ERs (∆) =

with
c b1 − 2
c b2 − 2
(b1 + b2 )1{0} (γ), m2 = −
(b1 + b2 )1{0} (γ),
1 − c 2b1
1 − c 2b2
1
c
(p−γ ∆2 )2 1[1,+∞) (γ) +
(c(b1 + b2 )1{0} (γ) + p−γ ∆2 1[0,1] (γ)),
v2 =
3
2(1 − c)
(1 − c)3

m1 =

e 2 by p−γ ∆2 .
where we approximate ∆
In the special case of n1 = n2 which leads to b1 = b2 = 2, we get


∆
ERs (∆) = Φ −hc
2



with
p
√
pmin(γ,1)/2−γ 1 − c p−γ ∆2
hc = p −γ 2 2
,
c(p ∆ ) 1[1,+∞) (γ)/2 + 4c1{0} (γ) + p−γ ∆2 1[0,1] (γ)
√
which is always smaller than one. Furthermore, for γ ∈ (0, 1) we get hc = 1 − c.
In Figure 2, we plot ERs (∆) as a function of ∆ ∈ [0, 100] for c ∈ {0.1, 0.5, 0.8, 0.95}.
19

0.0

0.1

0.2

0.3

0.4

0.5

ERp(∆)
ERs(∆), c = 0.1
ERs(∆), c = 0.5
ERs(∆), c = 0.8
ERs(∆), c = 0.95

0

20

40

60

80

100

Figure 2: Error rates ERp (∆) and ERs (∆) as functions of ∆ for c ∈ {0.1, 0.5, 0.8, 0.95}.
We also add the plot of ERp (∆) in order to compare the error rate of the two decision
rules. Since only finite values of ∆ are considered in the figure we put γ = 0 and also
1
1 +n2 −2
choose n1 = n2 . Finally, the ratio n1n+n
in the definition of a is approximated by 1−c
.
2 −p−1
We observe that ERs (∆) lies very close to ERp (∆) for c = 0.1. However, the difference
between two curves becomes considerable as c growths, especially for c = 0.95 and larger
values of ∆.

3.2

Finite-sample performance

In this section we present the results of the simulation study. The aim is to investigate
how good the asymptotic distribution of a linear combination of the discriminant function
coefficients θ̂ = lT â performs in the case of the finite dimension and of the finite sample
size. For that reason we compare the asymptotic distribution of the standardized θ̂ as
given in Theorem 4 to the corresponding exact distribution obtained as a kernel density
approximation with the Eppanechnikov kernel applied to the simulated data from the
standardized exact distribution which are generated following the stochastic representation of Corollary 1: (i) first, ξb , z0;b , ub are sampled independently from the corresponding

20

univariate distributions provided in Corollary 1; (ii) second, θ̂b is computed by using (11)
and standardized after that as in Theorem 4; (iii) finally, the previous two steps are repeated for b = 1, ..., B times to obtain a sample of size B. It is noted that B could be
large to ensure a good performance of the kernel density estimator.
In the simulation study, we take l = 1p (p-dimensional vector of ones). The elements
of µ1 and µ2 are drawn from the uniform distribution on [−1, 1] when γ > 0, while
the first ten elements of µ1 and the last ten elements of µ2 are generated from the
uniform distribution on [−1, 1] and the rest of the components are taken to be zero when
γ = 0. We also take Σ as a diagonal matrix, where every element is uniformly distributed
on (0, 1]. The results are compared for several values of c = {0.1, 0.5, 0.8, 0.95} and
the corresponding values of p, n1 , n2 . Simulated data consist of N = 105 independent
repetitions. In both cases γ = 0 and γ > 0 we plot two asymptotic density functions to
investigate how robust are the obtained results to the choice of γ.
In Figures 3-4, we present the results in the case of equal and large sample sizes (data
are drawn with γ = 0 in Figure 3 and with γ > 0 in Figure 4), while the plots in Figure
5 correspond to the case of one small sample and one large sample. We observe that the
impact of the incorrect specification of γ is not large, while some deviations are observed
in Figure 5 for small values of c. If c increases, then the difference between the two
asymptotic distributions becomes negligible. In contrast, larger differences between the
asymptotic distributions and the finite-sample one are observed for c = 0.8 and c = 0.95
in all figures, although their sizes are relatively small even in such extreme case.

References
Agarwal, A., Negahban, S., and Wainwright, M. J. (2012). Noisy matrix decomposition via
convex relaxation: Optimal rates in high dimensions. Annals of Statistics, 40(2):1171–
1197.
Bai, Z., Jiang, D., Yao, J.-F., and Zheng, S. (2009). Corrections to lrt on large-dimensional
covariance matrix by rmt. Annals of Statistics, 37(6B):3822–3840.
Bai, Z. and Silverstein, J. W. (2010). Spectral Analysis of Large Dimensional Random
Matrices. New York, NY: Springer Science+ Business Media, LLC.
Bickel, P. J. and Levina, E. (2004). Some theory for fisher’s linear discriminant function,’naive bayes’, and some alternatives when there are many more variables than
observations. Bernoulli, pages 989–1010.
Bodnar, T., Gupta, A., and Parolya, N. (2016). Direct shrinkage estimation of large
dimensional precision matrix. Journal of Multivariate Analysis, 146:223–236.

21

Bodnar, T., Gupta, A. K., and Parolya, N. (2014). On the strong convergence of the
optimal linear shrinkage estimator for large dimensional covariance matrix. Journal of
Multivariate Analysis, 132:215–228.
Bodnar, T., Mazur, S., and Okhrin, Y. (2017a). Bayesian estimation of the global minimum variance portfolio. European Journal of Operational Research, 256:292–307.
Bodnar, T., Mazur, S., and Parolya, N. (2017b). Central limit theorems for functionals of
large dimensional sample covariance matrix and mean vector in matrix-variate skewed
model. Scandinavian Journal of Statistics, under revision.
Bodnar, T. and Okhrin, Y. (2008). Properties of the singular, inverse and generalized
inverse partitioned wishart distributions. Journal of Multivariate Analysis, 99:2389–
2405.
Bodnar, T. and Okhrin, Y. (2011). On the product of inverse wishart and normal distributions with applications to discriminant analysis and portfolio theory. Scandinavian
Journal of Statistics, 38(2):311–331.
Bodnar, T. and Reiß, M. (2016). Exact and asymptotic tests on a factor model in low
and large dimensions with applications. Journal of Multivariate Analysis, 150:125–151.
Bodnar, T. and Schmid, W. (2008). A test for the weights of the global minimum variance
portfolio in an elliptical model. Metrika, 67(2):127–143.
Cai, T. and Liu, W. (2011a). Adaptive thresholding for sparse covariance matrix estimation. Journal of the American Statistical Association, 106(494):672–684.
Cai, T. and Liu, W. (2011b). A direct estimation approach to sparse linear discriminant
analysis. Journal of the American Statistical Association, 106(496):1566–1577.
Cai, T., Liu, W., and Luo, X. (2011). A constrained l1 minimization approach to
sparse precision matrix estimation. Journal of the American Statistical Association,
106(494):594–607.
Cai, T. T. and Jiang, T. (2011). Limiting laws of coherence of random matrices with
applications to testing covariance structure and construction of compressed sensing
matrices. Annals of Statistics, 39(3):1496–1525.
Chen, S. X., Zhang, L.-X., and Zhong, P.-S. (2010). Tests for high-dimensional covariance
matrices. Journal of the American Statistical Association, 105(490):810–819.
DasGupta, A. (2008). Asymptotic theory of statistics and probability. Springer Science
& Business Media.

22

Fan, J., Fan, Y., and Lv, J. (2008). High dimensional covariance matrix estimation using
a factor model. Journal of Econometrics, 147(1):186–197.
Fan, J., Liao, Y., and Mincheva, M. (2013). Large covariance estimation by thresholding
principal orthogonal complements. Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 75(4):603–680.
Fujikoshi, Y. and Seo, T. (1997). Asymptotic aproximations for epmcs of the linear and
the quadratic discriminant functions when the sample sizes and the dimension are large.
Random operators and stochastic equations, University of Toronto.
Givens, G. H. and Hoeting, J. A. (2012). Computational statistics. John Wiley & Sons.
Gupta, A. and Nagar, D. (2000). Matrix Variate Distributions. Chapman and Hall/CRC,
Boca Raton.
Gupta, A., Varga, T., and Bodnar, T. (2013). Elliptically contoured models in statistics
and portfolio theory. Springer, second edition.
Gupta, A. K. and Bodnar, T. (2014). An exact test about the covariance matrix. Journal
of Multivariate Analysis, 125:176–189.
Jiang, T. and Yang, F. (2013). Central limit theorems for classical likelihood ratio tests
for high-dimensional normal distributions. Annals of Statistics, 41(4):2029–2074.
Johnson, R. A., Wichern, D. W., et al. (2007). Applied multivariate statistical analysis.
Prentice hall Upper Saddle River, NJ.
Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components analysis. Annals of Statistics, 29(2):295–327.
Ledoit, O. and Wolf, M. (2003). Improved estimation of the covariance matrix of
stock returns with an application to portfolio selection. Journal of Empirical Finance,
10(5):603–621.
Mathai, A. and Provost, S. B. (1992). Quadratic Forms in Random Variables. Marcel
Dekker.
Muirhead, R. J. (1982). Aspects of Multivariate Statistical Theory. Wiley, New York.
Narsky, I. and Porter, F. C. (2013). Linear and quadratic discriminant analysis, logistic regression, and partial least squares regression. Statistical Analysis Techniques in
Particle Physics: Fits, Density Estimation and Supervised Learning, pages 221–249.
Rencher, A. C. (1998). Multivariate statistical inference and applications, volume 338.
Wiley-Interscience.
23

Rencher, A. C. and Christensen, W. F. (2012). Methods of multivariate analysis. John
Wiley & Sons.
Shao, J., Wang, Y., Deng, X., Wang, S., et al. (2011). Sparse linear discriminant analysis
by thresholding for high dimensional data. The Annals of statistics, 39(2):1241–1265.
Srivastava, M. S. and Kubokawa, T. (2007). Comparison of discrimination methods for
high dimensional data. Journal of the Japan Statistical Society, 37:123–134.
Tamatani, M. (2015). Asymptotic theory for discriminant analysis in high dimension low
sample size. Memoirs of the Graduate School of Science and Engineering, Shimane
University. Series B, Mathematics, pages 15–26.
Wyman, F. J., Young, D. M., and Turner, D. W. (1990). A comparison of asymptotic
error rate expansions for the sample linear discriminant function. Pattern Recognition,
23:775–783.

24

0.5

0.5

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

0

2

4

0.5

b) p = 250, n1 = n2 = 250

0.5

a) p = 50, n1 = n2 = 250

−2

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

c) p = 400, n1 = n2 = 250

−2

0

2

d) p = 475, n1 = n2 = 250

Figure 3: The kernel density estimator of the asymptotic distribution and standard normal
for θ̂ as given in Theorem 4 for γ = 0 and c = {0.1, 0.5, 0.8, 0.95}.

25

4

0.5

0.5

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

0

2

4

0.5

b) p = 250, n1 = n2 = 250

0.5

a) p = 50, n1 = n2 = 250

−2

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

c) p = 400, n1 = n2 = 250

−2

0

2

d) p = 475, n1 = n2 = 250

Figure 4: The kernel density estimator of the asymptotic distribution and standard normal
for θ̂ as given in Theorem 4 for γ > 0 and c = {0.1, 0.5, 0.8, 0.95}.

26

4

0.5

0.5

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

0

2

4

b) p = 250, n1 = 25, n2 = 475

0.5

0.5

a) p = 50, n1 = 25, n2 = 475

−2

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

0.0

0.0

0.1

0.2

0.3

0.4

Asymptotic True
Asymptotic False
Finite Sample

−4

−2

0

2

4

−4

c) p = 400, n1 = 25, n2 = 475

−2

0

2

d) p = 475, n1 = 25, n2 = 475

Figure 5: The kernel density estimator of the asymptotic distribution and standard normal
for θ̂ as given in Theorem 4 for γ > 0 and c = {0.1, 0.5, 0.8, 0.95}.

27

4

