DISTRIBUTIONS AND STATISTICAL POWER OF
OPTIMAL SIGNAL-DETECTION METHODS IN FINITE
CASES

arXiv:1702.07082v1 [] 23 Feb 2017

By Hong Zhang∗,‡ , Jiashun Jin§ and Zheyang Wu∗,†,‡
Worcester Polytechnic Institute‡ and Carnegie Mellon University§
In big data analysis for detecting rare and weak signals among n
features, some grouping-test methods such as Higher Criticism test
(HC), Berk-Jones test (B-J), and φ-divergence test share the similar
asymptotical optimality when n → ∞. However, in practical data
analysis n is frequently small and moderately large at most. In order
to properly apply these optimal tests and wisely choose them for
practical studies, it is important to know how to get the p-values
and statistical power of them. To address this problem in an even
broader context, this paper provides analytical solutions for a general
family of goodness-of-fit (GOF) tests, which covers these optimal
tests. For any given i.i.d. and continuous distributions of the input
test statistics of the n features, both p-value and statistical power
of such a GOF test can be calculated. By calculation we compared
the finite-sample performances of asymptotically optimal tests under
the normal mixture alternative. Results show that HC is the best
choice when signals are rare, while B-J is more robust over various
signal patterns. In the application to a real genome-wide association
study, results illustrate that the p-value calculation works well, and
the optimal tests have potentials for detecting novel disease genes
with weak genetic effects. The calculations have been implemented
in an R package SetTest and published on the CRAN.

1. Introduction. In big data analysis, signals are often buried within
a large amount of noises and are thus relatively weak and rare. It is ideal
to apply the optimal tests that are capable of detecting the minimal signals required by statistics. Through asymptotics, many theoretical studies
have brought exciting results on designing such optimal tests. In particular,
under the Asymptotically Rare and Weak (ARW) setting, the Higher Criticism (HC) and its various modifications, the Berk-Jones (B-J) type tests,
a spectrum of φ-divergence type tests, etc. were studied and proven asymptotically optimal [1, 2, 3, 4]. These optimal tests are attractive in many
scientific researches. For example, in large-scale genetic association studies,
a main strategy to find disease-associated genes is to determine whether
∗
†

Partially supported by NSF grant DMS-1309960.
Corresponding author.

1

2

H. ZHANG ET AL.

some of the genetic variants within candidate genes could affect disease outcome. Such genetic effects are often weak and rare, especially relative to the
cumulated noise level in big data.
However, for practical applications under finite cases the questions remain
on 1) how to analytically calculate p-values as well as statistical power, and
2) what are the real performances of those methods that are asymptotically
equivalent. First, to obtain the p-value for error control, the Monte-Carlo
simulation or permutation tests have significant limitations: (A) they require
daunting computation, and (B) empirical p-values are discrete, causing ties
among candidates. These issues are especially concerned when very small
p-values are demanded to handle a huge number of simultaneous tests. Secondly, for those optimal tests with the same asymptotic property under
n → ∞, it is important to understand their relative performance for finite
n. In order to solve these problems we need to calculate the distributions
of these tests under both H0 and H1 at each given n. Comparing with the
literature, this paper gives a complete answer by providing a comprehensive
calculation for a broad family of relevant tests.
1.1. Limitations of current methods. In general there are two types of
methods for distribution calculation. The first is to calculate the exact distribution. Recursive methods (e.g., Noe’s recursion[5, 6], Bolshev’s recursion
[7, 8], Steck’s recursion [9, 10], Ruben’s recursion [11], etc.) were developed to calculate the distribution of Kolmogorov-Smirnov type statistics
under H0 . In a similar fashion, Barnett and Lin [12] provided a calculation
method specifically for HC. Such recursive methods have heavy computation load, with complexity of O(n3 ). Moscovich, Nadler and Spiegelman [13]
reduced the computation to O(n2 ). However, all these methods assume the
full domain R = {1 ≤ i ≤ n} in the supremum formula of the relevant
tests (see the HC statistic in (3) as an example). However, calculation that
allows arbitrary R is important, because many tests improve performance
by properly restricting R. For example, a modified version of HC is under
R = {1 < i ≤ n/2, p(i) ≥ 1/n}. This is because that including too small
p-values p(i) could make HC heavy tailed; while including too big p-values
may unnecessarily reduces computational efficiency (c.f. [1]). Moreover, these
methods did not give statistical power calculation under H1 yet.
The second type methods are to approximate the distribution. KolmogorovSmirnov type statistics have been proven to converge in law to an extremevalue distribution [15, 16]. However, such convergence is too slow to be
accurate for even moderately large n [12]. Recently, Li and Siegmund (LS)
[14] developed an asymptotic approximation method. One problems is that

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES

3

Fig 1: Comparison among different methods for calculating the right-side
probability of the modified HC test (MHC) in (3) with R = {1 < i ≤
i.i.d.

n/2, p(i) ≥ 1/n}[1] under H0 : Xi ∼ N (0, 1). Black solid curves: by simulation; red dashed curves: by Corollary 1; green dotted curves: by Li and
Siegmund’s method[14].

this method cannot approximate the whole distribution. Figure 1 shows
that the LS calculation (green dotted curves) fails to match the simulation
(black solid curves) at small thresholds, whereas our calculation (e.g., the
red dashed curves) reveal the whole distributions. Moreover, comparing with
our study setting, LS assumes slightly more restriction for the supremum
searching domain R, and covers a narrower range of test types. Again, literature studies did not yet provide satisfiable method for deriving the relevant
distributions under H1 for statistical power calculation.
1.2. Our contribution. This paper has two folds of contributions. First,
it provides techniques for calculating the distributions of a general family
of goodness-of-fit (GOF) tests, which covers the optimal tests described
above. We give calculation methods for the exact as well as the approximated

4

H. ZHANG ET AL.

Fig 2: The association p-values for genes by exact calculation of four tests.
First row: HC 2004 and Berk-Jones; second row: reverse Berk-Jones and
HC 2008 .

distributions, balancing between accuracy and computation burden. The
methods allow (A) arbitrary truncation strategies (e.g., R for HC in (3)),
and (B) arbitrary null and alternative hypotheses as long as they are i.i.d.
and continuous. With such techniques, p-values and statistical power can be
calculated at any testing threshold.
Second, based on analytical calculation we systematically compared the
power of the asymptotically optimal tests. Simulations were applied for examing these tests under various n values and signals patterns in practice. We
also demonstrated the application of these weak-signal-sensitive tests in a
real genome-wide association study (GWAS) for detecting genes associated
to the Crohn’s disease (see Figure 2).

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES

5

In summary, the innovative contribution of this work lies on providing
a more mature statistical framework that revolves the problem of calculating distributions for a broader family of GOF tests under both H0 and
H1 . Thus, those asymptotically equivalent tests for detecting weak signals
can be apply, as well as compared by power, in practical applications. The
paper is organized as follows. In section 2 we formulate the problem under
a general setup of hypotheses and GOF tests, and brief the essence of our
strategy. The analytical results are presented in Section 3 for both exact
and approximated calculations. Through simulations Section 4 numerically
evidences the calculation accuracy, and provides power comparisons among
the asymptotically optimal tests. We show the application of the GOF tests
in a real GWAS in Section 5. In Section 6 we discuss relevant theoretical
and practical issues. The techniques of relevant proofs and lemmas are given
in the Appendix.
2. Problem Formulation and General Strategy.
2.1. Background on Detection of Weak and Sparse Signals. We consider
the signal detection problem through testing whether there exist “signals”
in noisy background. Under the broadly applicable Gaussian Means Model,
signals refer to the nonzero means in the normally distributed data values (or
statistics) X1 , X2 , ..., Xn . A test is said of detecting signals if it successfully
concludes the alternative when it is true. A typical setting for the null and
the alternative is
(1)

i.i.d.

H0 : Xi ∼ N (0, 1),

i.i.d.

H1 : Xi ∼ N (µ, 1) + (1 − )N (0, 1),

where  ∈ (0, 1) denotes the expected proportion of signals.
Under the Asymptotic Rare and Weak (ARW) setting, the
p proportion of
−α
signals is n = n , α ∈ (1/2, 1), and the mean is µn = 2r log(n), r ∈
(0, 1). Landmark studies [1, 17, 18] have provided the asymptotic detection
boundary in terms of a function curve of the signal strength and sparsity:

α − 1/2
1/2 < α ≤ 3/4
√
(2)
r = ρ? (α) =
3/4 < α < 1.
(1 − 1 − α)2
When the signals are below the curve, all tests will fail to distinguish H0
and H1 as n → ∞. Whenever signals are above the curve, the so-called
asymptotically optimal tests are able to make both the type I and the type II
error rates converge to zero. A particular optimal test is the Higher Criticism
(HC) test [19, 1]:
(3)

i/n − p(i)
√
HCn,R = sup n q
,
R
p(i) (1 − p(i) )

6

H. ZHANG ET AL.

where p(1) ≤ ... ≤ p(n) are the order statistics of pi = 2(1 − F0 (|Xi |)), i =
1, ..., n, the two-sided p-values of Xi . The supremum domain R regards to
the p-value magnitude, or the index i, or both. Note that in literature HC
formula could also be written as (e.g., [20, 12] )
P
i > t} − 2nΦ̄(t)
i {X
p
(4)
HC = sup
.
∗
t∈R
2nΦ̄(t)Φ(t)
Because the formula is monotone, the supremum domain R on p(i) is equivalent to the supremum domain R∗ on t.
Following that, a variety of modified Higher Criticism (HC), the BerkJones (B-J) type tests, a spectrum of φ-divergence type tests, etc. were
studied and proven asymptotically optimal [1, 2, 3, 4]. In this paper, we
provide calculation for the p-value and statistical power for relevant tests
under finite n.
2.2. Study Formulation. In this paper, we consider arbitrary i.i.d. and
continuous null and alternative hypotheses:
(5)

i.i.d.

H0 : Xi ∼ F0 ,

i.i.d.

Ha : Xi ∼ F1 , i = 1, ..., n.

We consider the goodness-of-fit (GOF) tests that evidence the distinction
between the data and any given null distribution. Consider either one-sided
p-values pi = 1 − F0 (Xi ) or two-sided p-values pi = 2(1 − F0 (|Xi |)). When
F0 is continuous, the null hypothesis is equivalent to
i.i.d.

H0 : pi ∼ Uniform(0, 1), i = 1, ..., n.
The typical idea of GOF is that if any Xi do not follow F0 , its p(i) shall
distinct from it expectation, which is roughly i/n. Following this idea we
consider a general family of tests,
(6)

i
Sn,R = sup f ( , p(i) ),
n
R

where the supremum domain
(7)

R = {i : k0 ≤ i ≤ k1 }

\

{p(i) : α0 ≤ p(i) ≤ α1 }

for given k0 ≤ k1 ∈ {1, ..., n} and α0 ≤ α1 ∈ [0, 1]. We assume that for fixed
x = i/n the function f (x, y) is monotonically decreasing in y = p(i) , so that
the smaller the input p-values, the larger the statistic and the stronger the
evidence against H0 .

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES

7

This GOF family contains lots of test statistics widely applicable to practice. For example, the simple one-sided Kolmogorov-Smirnov test statistic
(c.f. [21], page 447, denoted KS + here) is a classic GOF test, which directly measures the difference between p(i) and i/n. That is, the f function
is defined by
fKS + (x, y) = x − y.

(8)

To improve the test, the difference between p(i) and i/n should be scaled
with regard to p(i) or i/n. This is because smaller p(i) values are more important to evidence agains H0 . Such scaled KS tests are related to the Higher
Criticism (HC) statistics proposed in 2004 and 2008 [1, 2], respectively,
where the f functions are defined as
√

x−y
np
;
y(1 − y)
√
x−y
.
fHC 2008 (x, y) = n p
x(1 − x)

fHC 2004 (x, y) =
(9)

Jager and Wellner introduced a collection of φ-divergence statistics [4],
which are also based on the supremum of functions:
1
(1 − xs y 1−s − (1 − x)s (1 − y)1−s ), s 6= 0, 1,
s(1 − s)
x
1−x
f1φ (x, y) = x log( ) + (1 − x) log(
),
y
1−y
y
1−y
f0φ (x, y) = y log( ) + (1 − y) log(
).
x
1−x
fsφ (x, y) =
(10)

For certain s (e.g., s = 2 or −1) these statistics are two-sided in the sense
that switching the values of x = i/n and y = p(i) gives the same statistic.
However, this property is not appropriate in the scenario of signal detection.
It is the relationship p(i) < i/n, instead of the opposite, indicates signals.
Thus for signal detection purpose, one-sided test is more reasonable and
more powerful, as is the same idea for defining KS and HC. Thus, here
we consider the one-sided version of φ-divergence statistics, which can be
achieved by a simple adjustment of the f function to be, for example,
 q
 2nf φ (x, y)
y ≤ x,
q s
(11)
fs (x, y) =
 − 2nf φ (x, y)
y > x.
s

Now for all s, fs (x, y) is guaranteed decreasing in y. Such one-sided φdivergence statistics cover HC exactly: f2 = fHC 2004 and f−1 = fHC 2008 .

8

H. ZHANG ET AL.

Also, s = 1 and 0 correspond to the Berk-Jones statistic [22, 1, 14] and the
reverse Berk-Jones statistic [22], respectively. Note that the input p-values
themselves could be two-sided in order to accommodate the consideration
of signal directionality in practical problems.
2.3. Calculation Strategy. In the following we first introduce the essential
idea of the calculation. Under various settings and assumptions, detailed
strategies for obtaining the exact and approximated distributions will be
described in Section 3.
For any given continuous CDFs in (5), we define a monotone transformation function in (0, 1):

x
under H0 : F0 ,
(12)
D(x) =
under H1 : F1 6= F0 .
1 − F1 (F0−1 (1 − x))
Note that for any p-value pi , we have D(pi ) ∼ U nif orm(0, 1) under either
H0 or H1 .
Consider the function f (x, y) of any statistic Sn,R in (6), for each fixed x
define its inverse function
(13)

g(x, ·) = f −1 (x, ·).

For example, for the HC statistics defined (9), the g functions are
√ √
b2 /n−(b/ n) b2 /n+4x(1−x)
1
gHC 2004 (x, b) = 1+b2 /n [x +
];
2
(14)
√ p
gHC 2008 (x, b) = x − (b/ n) x(1 − x).
In general if the closed form of g(x, ·) is not available, it can always be found
numerically, since f (x, y) is strictly decreasing in y.
Now under either H0 or H1 , the CDF function of Sn,R is

(15)

i
P (Sn ≤ b) = P (sup f ( , p(i) ) ≤ b)
n
R
\
i
= P ( {p(i) > g( , b)})
n
R

= P {D(p(i) ) > D(g( ni , b)), all i and p(i) in R}.
The key idea is that under either H0 or H1 , U(i) := D(p(i) ) is the ith order
statistic of Uniform(0, 1), and the joint distribution of U(i) , i = k0 , ...., k1 ,
can be studies one way or another for getting the final probability.
To simplify the presentation, we list below the notations to be referred
later on.

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES

9

(N1) uk := D(g( nk , b) ∨ α0 ), following the definitions in (12)–(15), and a
potential constant α0 ≥ 0 in (7).
(N2) F̄B(α,β) (x) denotes the survival function of Beta(α, β) distribution.
(N3) FΓ(α) (x) and F̄Γ(α) (x) denote the CDF and survival function of Gamma(α, 1)
distribution, respectively, where the shape parameter is α, the scale
parameter is 1.
(N4) Based on the notation (N3) define
hk (x) := xFΓ(k−1) (kx) − FΓ(k) (kx).
(N5) fP (λ) (x) denotes the probability mass function of P oisson(λ) distribution.
3. Analytical Results.
3.1. Exact Calculations. The first theorem provides the exact calculation
for the distribution of GOF in (6), where the supremum domain R involves
truncation of the index i. For example, the initial HC was defined with
R = {1 ≤ i ≤ n/2} [1].
Theorem 3.1. Consider any GOF statistic in (6) with R = {k0 ≤ i ≤
k1 } for given 1 ≤ k0 ≤ k1 ≤ n. Let m = n − k1 + 1. Following the notations
(N1), (N2), and
ak1 =
ak =

n!
(n−k1 +1)! F̄B(1,m) (uk1 ),

and

kX
1 −k
ujk+j−1
n!
F̄B(k1 −k+1,m) (uk1 ) −
ak+j ,
(n − k + 1)!
j!

k = k1 − 1, ..., 1.

j=1

Under either H0 or H1 we have
P (Sn ≤ b) = F̄B(k1 ,m) (uk1 ) −

kX
1 −1
i=k0

uii
ai+1 .
i!

Another type of truncation is based on the value of p(i) . In the case of HC
under the null of N (0, 1), this truncation is equivalent to the truncation on
t in (4). The following theorem gives the exact calculation for the general
GOF tests in (6) with such truncations.
Theorem 3.2. Consider any GOF with statistic in (6) with R = {α0 ≤
p(i) ≤ α1 } for given 0 ≤ α0 < α1 ≤ 1. Following the notations (N1) and

10

H. ZHANG ET AL.

(N2), define
β0 = D(α0 ),
aj (k) =

β1 = D(α1 ),

β0i−1 (1−β1 )n−j+1
(i−1)!(n−j+1)!
j−k l
X
uk+l−1
uj−1

cij =

n!
)−
β1j−k F̄B(j−k,1) (
(j − k)!
β1

aj (j) = 0,

l!

l=1

aj (k + l), and

1 ≤ i ≤ k1 , i < j ≤ n + 1, k = 1, ..., j − 1

Under either H0 or H1 , we have
P (Sn,R ≤ b) =

k1 n+1
X
X

cij aj (i)

i=1 j=i+1

The most general R is in (7), which defines the truncation for both the
index and the p-values. The following theorem provides the calculation for
exact distribution under this general setup, for which no literature has provided solution before.
Theorem 3.3. Consider any GOF with statistic in (6)with R = {α0 ≤
p(i) ≤ α1 } ∩ {k0 ≤ i ≤ k1 } for given 1 ≤ k0 ≤ k1 ≤ n and 0 ≤ α0 < α1 ≤ 1.
Following the notations in Theorem 3.2 and (N1) and (N2), define
(
(
i
i ≥ k0
j
j ≤ k1 + 1
ĩ =
, j̃ =
, β̃0 = β0 I{i<k0 } ,
k0 i < k0
k1 + 1 j > k1 + 1
j̃−k

aj (k) =

X ulk+l−1
uj̃−1
n!
(j−k)
β1
)−
aj (k + l), and
F̄B(j̃−k,j−j̃+1) (
(j − k)!
β1
l!
l=1

aj (j̃) = 0,

1 ≤ i ≤ k1 , ĩ < j ≤ n + 1, k = 1, ..., j̃ − 1.

Under either H0 or H1 , we have
P (Sn,R ≤ b)
=

k1 n+1
X
X
i=1 j=ĩ+1




j̃−1
j−i
k−i+1
X
˜
u
−
β̃
0
n!(β
−
β̃
)
(u
−
β
)
j̃−1
1
0
0
k
cij 
)−
aj (k + 1) .
F̄B(j̃−i,j−j̃+1) (
(j − i)!
(k
−
i
+ 1)!
β1 − β̃0
k=ĩ

The study setting in Li and Siegmund [14] is a special case of the general
truncation with α1 = 1. For example, the “modified HC” is to take R =
{1 < i ≤ n/2, p(i) ≥ 1/n}[1]. Li and Siegmund’s calculation provided an
approximation for the modified HC under H0 . As shown in Figure 1, the
approximation does not give the whole distribution. For this setting, the
corollary below gives the exact distribution for the general GOF in (6) under
both H0 and H1 .

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 11

Corollary 1. Consider any GOF with statistic in (6) with R = {α0 ≤
p(i) } ∩ {k0 ≤ i ≤ k1 } for given 1 ≤ k0 ≤ k1 ≤ n and α0 > 0. Following the
notations in Theorem 3.1, 3.2 and (N1) and (N2), define
ci =

β0i−1
,
(i − 1)!

1 ≤ i ≤ k1 .

Under either H0 or H1 , we have


k1
kX
1 −1
n+1−i
k+1−i
X
n!(1
−
β̃
)
u
−
β̃
(u
−
β̃
)
0
0
0
k
k
P (Sn,R ≤ b) =
ci 
F̄B(k1 +1−i,m) ( 1
)−
ak+1  .
(n + 1 − i)!
(k + 1 − i)!
1 − β̃0
i=1

k=ĩ

To sum up, Theorem 3.3 covers Theorem 3.1 by fixing i = 1, j = n + 1
(so that cij = 0) and α0 = 0, α1 = 1. It covers Theorem 3.2 by letting
k0 = 1, k1 = n, and covers Corollary 1 by fixing j = n + 1 and α1 = 1.
Regarding the computational complexity, the calculation given by Theorem 3.1 is equivalent to solving an upper triangular linear systems using
backward substitution, which is O(n2 ), an improvement from O(n3 ) required
by other exact distribution calculations [5, 6, 12]. More importantly, our
method can handle a more general family of GOFs defined in (6) with a
more flexible supremum domain R. When R is more complicated, Corollary
1 is still O(n2 ) because the inner loop ak is not dependent on i. Theorem
3.2 and Theorem 3.3 are O(n3 ).
3.2. Approximate the distributions. In this section, we provide several
calculation methods to approximate the distributions of the GOF family.
Under more restricted conditions, computation load could be significantly
reduced. For simplicity we provide results here for the supremum domain
R = {k0 ≤ i ≤ k1 }. More general results for R in (7) can be obtained by
following the similar idea of the calculation in Theorem 3.3.
The theorem below gives an approximated calculation based on a joint
Gamma distribution.
Theorem 3.4. Consider any GOF with statistic in (6) with R = {k0 ≤
i ≤ k1 }. Following the notations (N1) and (N3), define
dk = (n + 1)D(g( nk , b)), k = k0 , ..., k1 .
k−1 j
X
dk1 −k+j
ck = F̄Γ(k) (dk1 ) −
ck−j , k = 2, ..., k1 , and
j!
j=1

c1 = F̄Γ(1) (dk1 ).

12

H. ZHANG ET AL.

Under either H0 or H1 , we have
P (Sn ≤ b) = (1 + o(1))(F̄Γ(k1 ) (dk1 ) −

kX
1 −1
k=k0

dkk
ck −k ).
d! 1

To further reduce the computation, we can also deduce a one-step formula
under stronger assumptions. In particular, if D(g( nk , b)) is a linear function
of k, we can provide a closed-form formula for the distribution that gives
the same accuracy as the above theorem.
Theorem 3.5. Consider any GOF with statistic in (6) with R = {1 ≤
i ≤ k1 }. Assume the function D(g( nk , b)) = a+λk, for some λ ≥ 0. Following
the notation (N3) and (N4), under either H0 or H1 we have
P (Sn ≤ b) = (1 + o(1))e−a (1 − λ + hk1 (λ)).
One example that satisfies the linearity of D(g( nk , b)) is the simple KolmogorovSmirnov (KS + ) test in (8) under H0 , where a = −(n + 1)b and λ = n+1
n . In
this case, based on the above theorem we have the following corollary
Corollary 2. Consider the test statistic KS + in (8) under H0 . Following the notation (N3) and (N4), for b ≤ n1 , we have
P (KS + ≤ b) = (1 + o(1))e(n+1)b (−

1
n+1
+ hk1 (
)).
n
n

The requirement of linear D(g( nk , b)) in Theorem 3.5 is stringent. However,
if D(g( nk , b)) is close to linear, we can still further simplify the calculation
of Theorem 3.4. As one example, Theorem 3.6 below provides such a calculation formula for GOF and the conditions on D(g( nk , b)). The strategy for
proving this theorem largely follows the idea in Li and Siegmund [14]. The
key difference is that, instead of using the beta distribution in Li and Siegmund, we use the gamma distribution, which has a simpler density function
for addressing a wider family of the GOF tests. See Appendix for details.
Theorem 3.6. Consider GOF tests in (6) with R = {k0 ≤ i ≤ k1 }.
Follow notations (N1) – (N5), and define dk = (n + 1)D(g( nk , b)), d0k =
√
d
(n + 1) dx
D(g( nk , b)), and k ∗ = min{k1 − k, n}. Assume D(g(x, b)) satisfies
1. D(g(x, b)) < 1 is increasing and convex in x for
d
2. dx
D(g(x, b)) < 1, and

k0
n

≤x≤

k1
n,

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 13

3. D(g(k/n, b)) <

k
n+1 ,

for k > 1 and large n.

Under either H0 or H1 we have
P (Sn ≥ b) = (1 + o(1))

k1
X
k=k0

(1 −

d0k
d0
+ hk∗ ( k ))fP (dk ) (k).
n
n

Such a close-to-linear property of D(g( nk , b)) can be satisfied by HC 2004
under H0 . In this case, D(g( nk , b)) = g(x, b) is given in (14), and the condi√
tions can be satisfied when b is in the order of O( n). Thus, we can get a
computation-easy formula.
Corollary 3.
Consider the test statistic HC 2004 in (9) with R =
{k0 ≤ i ≤ k1 }. Let b0 = √bn be a positive constant > 2x − 1, kn0 < x < kn1 .
Define function
q
1
2
2
g(x, b0 ) =
[x
+
(b
−
b
0 b0 + 4x(1 − x))/2],
0
1 + b20
1
b0 (1 − 2x)
g 0 (x, b0 ) =
[1 − p 2
].
2
1 + b0
b0 + 4x(1 − x)
Following the notation (N2), under H0 , we have

k1 
X
0 k
2004
0 k
P (HC
≥ b) = (1+o(1))
1 − g ( , b0 ) + hk∗ (g ( , b0 )) fP (g( k ,b0 )n) (k).
n
n
n
k=k0

The above formula is different from that given in Li and Sigmund [14].
However, for the HC under H0 , both formulas require the threshold b =
√
O( n). Thus, both methods do not get the whole distribution. Furthermore, the accuracy depends on the linear approximation of the D(g( nk , b))
function, which could be far from the truth under general H1 . Thus this
type of calculation may not be very attractive for calculation statistical
power. Meanwhile, numeric results in Section 4 show that both methods are
accurate for calculating small p-values under relatively large threshold b.
4. Numerical results. In this section we first evidence the accuracy
of our methods by comparing the calculations with the Monte-Carlo simulations under various settings of H0 and H1 . Then, based on calculation we
compare the finite-n performance of the asymptotically optimal tests over
various signal patterns. The supremum domain is R = {1 ≤ i ≤ n/2} by default, and will be specified otherwise. The simulations run 5,000 repetitions
by default.

14

H. ZHANG ET AL.
i.i.d.

4.1. Evaluate the accuracy of calculations. Under H0 : Xi ∼ N (0, 1), i =
1, ..., n. Figure 3 shows the right-tail probability of HC over the threshold
b. Comparing with simulation (black solid curves), the exact calculation by
Theorem 3.1 (cyan dashed curves) has a perfect match. The approximation
by Theorem 3.4 is accurate over the whole distribution when n is fairly big.
The closed-formula calculation methods by Li and Siegmund [14] (blue dotted curves) and by Corollary 3 (green dashed curves) do not fit the whole
distribution curve, but both can provide good approximation for calculation
small p-values at large threshold.
For distribution calculation under the alternative, we considered the null
of i.i.d. N (0, 1) (on which the input p-values are based) and the alternatives:
i.i.d.

H1 : Xi ∼ (1 − )N (0, 1) + N (1, 1), or
i.i.d.

H1 : Xi ∼ (1 − )N (0, 1) + tν .
Figure 4 shows the right-tail probability of HC under the normal mixture
alternative with µ = 1,  = 0.1 (row 1), or under the Student’s t distribution with degrees of freedom ν = 5 (row 2). In both cases the distribution
curves by exact calculation (Theorem 3.1, cyan dashed curves) and by approximation (Theorem 3.4, red dot-dashed) are close to simulation (black
solid curves). As expected, the approximation is more accurate for larger n.
Besides the normal distributions, Theorems 3.1 – 3.4 can handle any given
continuous F0 and F1 . Here we test on four settings studied in the initial
paper of HC [1]. The first setting regards a Chi-squared model:
i.i.d.

H0 : Xi ∼ χ2ν (0),

vs.

i.i.d.

H1 : Xi ∼ (1 − )χ2ν (0) + χ2ν (δ).

where ν is the degree of freedom, δ is the non-centrality parameter. The
second setting is a Student’s t mixture model:
i.i.d.

H0 : Xi ∼ tν (0),

vs.

i.i.d.

H1 : Xi ∼ (1 − )tν (0) + tν (δ).

The third setting is a chi-squared-exponential mixture model,
i.i.d.

H0 : Xi ∼ exp(ν),

vs.

i.i.d.

H1 : Xi ∼ (1 − ) exp(ν) + χ2ν (δ).

The fourth setting concerns a generalized normal distribution (also known
as power exponential distribution) model,
i.i.d.

H0 : Xi ∼ GNp (0, σ),

vs.

i.i.d.

H1 : Xi ∼ (1 − )GNp (0, σ) + GNp (µ, σ),

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 15

Fig 3: Comparison among different calculations for the distributions of HC
i.i.d.
over various n value under H0 : Xi ∼ N (0, 1). Black solid curves: by
simulations; cyan dashed curves: by Theorem 3.1 for exact distribution; red
dot-dashed curves: by Theorem 3.4; blue dotted curves: by Li and Siegmund
[14]; green dashed curves: by Corollary 3 .

16

H. ZHANG ET AL.

i.i.d.

Fig 4: The alternative distributions of HC under H0 : Xi ∼ N (0, 1) vs.
H1 : Xi ∼ 0.9N (0, 1) + 0.1N (1, 1) (row 1), or H1 : Xi ∼ 0.5N (0, 1) + 0.5t5
(row 2). Column 1: n = 10; column 2: n = 100. Black solid curves: by
simulations; cyan dashed curves: by Theorem 3.1; red dot-dashed curves: by
Theorem 3.4.

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 17

where the probability density function of GNp (µ, σ) is
1
|x − µ|p
exp(−
),
Cp
pσ p

Cp = 2p1/p Γ(1 + 1/p)σ.

Notice that GN1 (µ, σ) is the Laplace distribution and GN2 (µ, σ) is N (µ, σ 2 ).
Each row of Figure 5 illustrates the alternative distribution of HC under each
of the four settings for n = 10 (left column) and 100 (right column). Again
the distribution calculation is fairly accurate in all cases.
For the closed-form calculation formula given by Theorem 3.5, the boundary is assumed linear: D(g( ni , b)) = a + λk ≥ 0 in (15). One example is the
KS + in (8) under H0 . Figure 6 demonstrates the accuracy of the calculation
based on either fixed slope λ = 0.5 or fixed intercept a = 0.5. Here k0 = 1,
k1 = n = 50. As the boundary a + λk increases, the probabilities from both
calculation and simulation decrease and well-matched as expected.
4.2. Compare statistical power of asymptotically optimal tests. All φdivergence statistics with s ∈ [−1, 2] possess the same asymptotic optimality property for detecting weak and sparse signals [4]. It is of interest
to know the performance of such statistics under finite n. Here we study
s = 2, 1, 0, −1, which correspond to HC 2004 , the Berk-Jones statistic, the
reverse Berk-Jones statistic, and HC 2008 , respectively.
To show the calculation accuracy, we calculated (by Theorem 3.1) the
critical values at the significance levels of 10%, 5% and 1%. Then at these
critical values we got the empirical type I error rates through simulation
(10,000 repetitions). Table 1 presents the thresholds by calculation and the
correspondingly empirical type I error rates. The closeness of these empirical
and nominal type I error rates evidences that the calculation for these tests
is accurate.
Now through the calculation based on Theorem 3.1, we systematically
compared the power of these tests under
i.i.d.

H0 : Xi ∼ N (0, 1),

vs.

i.i.d.

H1 : Xi ∼ (1 − )N (0, 1) + N (µ, 1).

With the type I error rate controlled at 5%, Figure 7 provides the statistical
power at various µ, n and . There are a few interesting observations. First,
HC 2004 performs well when signals are sparse. However, it seems more relevant to the average number of signals, i.e., n, rather than the proportion
 itself. For example, at fixed n = 5 (panels in the first column), HC 2004
is always the best, while at fixed  = 0.05 (panels in the diagonal), HC 2004
becomes worse when n increases. Second, Smaller s parameter (s = −1 for

18

H. ZHANG ET AL.

Fig 5: The alternative distributions of HC under four non-normal settings
for H0 and H1 . Column 1: n = 10; column 2: n = 100. Black solid curves:
by simulations; cyan dashed curves: by Theorem 3.1; red dot-dashed curves:
by Theorem 3.4.

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 19

Fig 6: Probability in (15) with boundary D(g( ni , b)) = a + λk. Black solid
curves: by simulations; red dot-dashed curves: by Theorem 3.5.

Table 1
Type I error rates: exact calculation vs. simulation.
s
2

1

0

-1

n
10
50
100
10
50
100
10
50
100
10
50
100

10%
threshold simulation

threshold

3.357
3.507
3.539
2.181
2.408
2.478
1.750
2.040
2.136
1.618
1.909
2.010

4.648
4.714
4.723
2.504
2.716
2.780
1.974
2.301
2.402
1.838
2.165
2.271

0.992
0.102
0.103
0.101
0.098
0.104
0.100
0.101
0.101
0.098
0.099
0.107

5%
simulation
0.049
0.050
0.049
0.050
0.048
0.049
0.049
0.047
0.051
0.051
0.049
0.052

threshold
10.088
10.102
10.102
3.110
3.300
3.354
2.390
2.803
2.915
2.227
2.662
2.777

1%
simulation
0.010
0.011
0.009
0.011
0.010
0.009
0.011
0.011
0.010
0.009
0.009
0.010

20

H. ZHANG ET AL.

HC 2008 and s = 0 for the reverse Berk-Jones) gives better power for denser
signals. However, when signals are sparse, they are less powerful because
these statistics are insensitive to small p-values as they are weighted by i/n
rather than p(i) . Third, Berk-Jones statistic has a more robust performance
over various µ, n and , which is consistent with the finding of Li and Siegmund [14].
Moreover, it is interesting to see the performance of these methods along
the asymptotic detection boundary. This is because the optimal tests are
most valuable for detecting subtle signals around the detection boundary,
for which sub-optimal tests will have asymptotically zero power. Here we
consider the ARW setting with the detection boundary given in (2). Figure
8 shows the power (calculated by Theorem 3.1) of the four methods over
the sparsity parameter α ∈ (1/2, 1) and r calculated accordingly. It shows
that the power curves of HC 2004 and B-J have the similar pattern, while
the performance of HC 2008 and reverse B-J are alike. In the less sparse case
with α ∈ (1/2, 3/4), B-J is preferred; in the sparser case with α ∈ (3/4, 1),
HC 2004 outperforms the others.
Last but not least, the supremum domain R in (6) is quite relevant in
constructing the test statistic. Here we compare HC 2004 under R = {1 ≤
i ≤ n/2} with the modified HC under R = {1 < i ≤ n/2, p(i) ≥ 1/n}. Figure
9 shows that the MHC performs poorly when the number of signals is small,
whereas it improves the performance when the number of signals increases.
When signals are sparse, MHC is less powerful because it tends to exclude
signals by considering only those p-values bigger than 1/n, the latter is a
fairly large value when n is small. Thus, in practice when n is not too big,
there may be no need to truncate p-values by 1/n.
5. A Genome-wide Association Study for Crohn’s Disease. According to the genetics of complex diseases, disease-associated markers usually have moderate to small effects[23]. In genome-wise studies that tend
to screen as many markers as possible, the number of true disease markers
often account for a small proportion of the total candidates. Therefore, it
is appealing to apply tests that are asymptotically optimal for detecting
rare and weak signals. In this section, based on a logistic regression model,
we applied optimal tests to a SNP-set association study at the gene level.
Specifically, for the ith individual, i = 1, ..., N , let yi = 1 (or 0) indicate
the case (or control), let Gi = (Gi1 , ..., Gin ) denote the genotype vector of
n SNPs in a given gene, and let Zi = (1, Zi1 , Zi2 ) contain the intercept
and the first two principal components of the genotype data, which serve
the purpose for controlling potential population structure [24]. The logis-

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 21

Fig 7: Comparison of statistical power (at type I error rate 5%).

22

H. ZHANG ET AL.

Fig 8: Statistical power along the ARW detection boundary (at type I error
rate 5%).

Fig 9: Power comparison for HC vs. MHC (at type I error rate 5%).

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 23

tic regression model is logit(E(yi |Gi , Zi )) = G0i β + Zi0 γ, where β and γ are
the coefficient parameters. A classic marginal-test statistic of the jth SNP,
adjusted for the non-genetic measures, is [25, 26]
Uj =

N
X

Gij (yi − ỹi ), j = 1, ..., n,

i=1

where ỹi is the fitted outcome value under H0 : β = 0, i.e., none SNPs in the
gene are associated. Under H0 , the vector of statistics U = (U1 , ..., Un )
N (0, Σ), where Σ can be estimated by
Σ̂ = G0 W G − G0 W Z(Z 0 W Z)−1 Z 0 W G,
where G = (Gij ) and Z = (Zij ) are the corresponding design matrices,
and W = diag(ỹi (1 − ỹi )). After de-correlation we get test statistics X =
1
Σ̂− 2 U
N (0, In×n ). Thus the i.i.d. condition of H0 in (5) is reasonably
satisfied, and our calculation methods can apply to obtain the p-value of a
GOF statistic in (6), which measures the association of the given gene.
We examined HC 2004 , B-J, reverse B-J, and HC 2008 (again, they corresponding to the φ-divergence statistics with s = 2, 1, 0, −1, respectively).
The GWAS data from NIDDK-IBDGC (National Institute of Diabetes, Digestive and Kidney Diseases - Inflammatory Bowel Disease Genetics Consortium) contain 1,145 individuals from non-Jewish population (572 Crohn’s
disease cases and 573 controls) [27]. After typical quality control for genotype data, 308,330 somatic SNPs were grouped into 15,857 genes according
to their physical locations.
Figure 2 gives the QQ plots of the gene-association-indicating p-values
calculated by Theorem 3.1. The genomic inflation factors, i.e., the ratios of
empirical median of -log(p-values) vs. the expected median under H0 , are
all close to 1, evidencing that the genome-wide type I errors were well controlled. Among the four tests, the B-J seemed having higher power because
it yielded more genes significantly above the red line of the H0 -expected
p-values. Among the top ranked genes, IL23R and CARD15 (also known as
NOD2) are well-known Crohn’s disease genes [28, 29, 27]. Gene NPTX2 was
top ranked by both HC 2004 and B-J. It hasn’t been reported to be directly
associated, but it encodes a member of the family of neuronal petraxins,
synaptic proteins that are related to C-reactive protein [30], while C-reactive
protein is an indicator for the activity level of Crohn’s disease [31]. Furthermore, NPTX2 has an important paralog gene APCS (www.genecards.org),
which is related to arthritis, a disease highly correlated with Crohn’s disease [32]. Gene SLC44A4 is also related to the pathophysiology of Crohn’s.

24

H. ZHANG ET AL.

Defects in this gene can cause sialidosis, a lysosomal storage disease [30]
that results from a deficiency of the enzyme sialidase, the later is important
for various cells to defend against infection [33]. Gene BMP2 was identified by B-J, reversed B-J, and HC 2008 . This gene could also be relevant
because it is associated with digestive phenotypes, especially colon cancer
[34, 35]. Meanwhile, for those top ranked genes, further studies are needed
to validate.
6. Discussion. This paper provides techniques to calculate the exact
and approximated null and alternative distributions of a general family of
GOF tests. Thus we can calculate both the p-value and the statistical power
of these tests. These calculations are important for both practical applications of GOF tests and for performance comparisons under finite n. Comparing with the relevant literature, our methods are novel, accurate, general
to broad statistic family and supremum domains.
To calculate the exact distribution, the result in Theorem 3.1 brings down
the computational complexity from O(n3 ) to O(n2 ) when comparing with
corresponding literature methods. In the meanwhile, when n and the search
range k1 −k0 are large, the calculation could suffer from the loss of significant
digits. In this case, we could truncate the summation to the first 25 - 30
terms, which yields a fairly accurate result and saves computation time too.
Moreover, we could also apply proper approximated calculations, for which
the accuracy increases together with n.
APPENDIX A: PROOFS
A.1. Proofs of the Main Theorems.
A.1.1. Proof of Theorem 3.1. For k = k1 , ..., 1, define
Z xk
Z xk+1
Z 1
1
n!
ak =
(1 − xk1 )n−k1
...
dxk ...dxk1 −1 dxk1 .
uk (n − k1 )!
uk −1
uk
1

1

n!
(n−k1 +1)! F̄B(1,m) (uk1 ),

Then obviously ak1 =
and for k ≤ k1 − 1,
Z 1
Z xk
Z xk+2
1
n!
n−k1
ak =
(1 − xk1 )
...
xk+1 dxk+1 ...dxk1 −1 dxk1 − uk ak+1
uk1 (n − k1 )!
uk1 −1
uk+1
Z 1
kX
1 −k
ujk+j−1
xkk11 −k
n!
n−k1
=
(1 − xk1 )
dxk1 −
ak+j
(k1 − k)!
j!
uk (n − k1 )!
j=1

1

=

n!
F̄
(uk ) −
(n − k1 + k)! B(k1 −k+1,m) 1

kX
1 −k
ujk+j−1
j=1

j!

ak+j

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 25

Now by Lemma 1,
P (Sn ≤ b) = P {U(k) > uk , k0 ≤ k ≤ k1 }
Z 1
Z xk +1
Z xk
xkk00 −1
0
1
n!
n−k1
=
...
(1 − xk1 )
dxk0 ...dxk1 −1 dxk1
(k0 − 1)!
uk0
uk1 (n − k1 )!
uk1 −1
Z xk +2 k0
Z xk
Z 1
xk0 +1
uk0
0
1
n!
n−k1
...
(1 − xk1 )
dxk0 +1 ...dxk1 −1 dxk1 − k0 ak0 +1
=
k0 !
k0 !
uk0 +1
uk1 −1
uk1 (n − k1 )!
Z 1
kX
1 −1
xk1 −1
n!
uii
=
(1 − xk1 )n−k1 k1
dxk1 −
ai+1
(k1 − 1)!
i!
uk (n − k1 )!
i=k0

1

kX
1 −1

=F̄B(k1 ,m) (uk1 ) −

i=k0

uii
ai+1 .
i!

A.1.2. Proof of Theorem 3.2. The idea is to use total probability theorem.
P (Sn,R ≤ b) =

k1 n+1
X
X


P {Sn,R ≤ b} ∩ {exactly p(i) ...p(j−1) fall in [α0 , α1 ]}

i=1 j=i+1

=

k1 n+1
X
X

Pij

i=1 j=i+1

Notice that 1 ≤ i ≤ n, i + 1 ≤ j ≤ n + 1, we have
Pij = P (U(i−1) < β0 , U(i) ≥ ui , ..., U(j−1) ≥ uj−1 , U(j−1) ≤ β1 , U(j) > β1 ),
where the joint density of U(i−1) , ..., U(j) is
f (xi−1 , ..., xj ) =

n!
xi−2 (1 − xj )n−j ,
(i − 2)!(n − j)! i−1

0 ≤ xi−1 ≤ ... ≤ xj ≤ 1.

Then
i−2
xi−1
Pij = n!
···
dxi−1 dxi ...dxj−1 dxj
(i − 2)!
β1
uj−1
ui
0
Z 1
Z β1
Z xi+1
Z β0
i−2
xi−1
(1 − xj )n−j
= n!
dxj
···
dxi ...dxj−1
dxi−1
(n − j)!
(i − 2)!
β1
uj−1
ui
0
Z β1
Z xi+1
= cij n!
···
dxi ...dxj−1 .

Z

1

(1 − xj )n−j
(n − j)!

uj−1

Z

β1

Z

xi+1

Z

β0

ui

Direct calculation similar to the proof of Theorem 3.1 gives the final result.

26

H. ZHANG ET AL.

A.1.3. Proof of Theorem 3.3. Following the idea in the proof of Theorem
3.2,
P (Sn,R ≤ b) =

k1 n+1
X
X


P {Sn,R ≤ b} ∩ {exactly p(i) ...p(j−1) fall in [α0 , α1 ]}

i=1 j=ĩ+1

=

k1 n+1
X
X

Pij .

i=1 j=ĩ+1

For each feasible pair of (i, j), direct calculation shows Pij can be concisely
written as
(β1 − xj̃−1 )j−j̃

β1

Z
Pij = cij n!

(j − j̃)!

uj̃−1


= cij n!

Z

)j−i

(β1 − β̃0
(j − i)!

xj̃−1

Z
···

uj̃−2

xĩ+1

(xĩ − β̃0 )ĩ−i
dxĩ ...dxj̃−2 dxj̃−1
(ĩ − i)!

uĩ

F̄B(j̃−i,j−j̃+1) (

uj̃−1 − β̃0
β1 − β̃0

)−

j̃−1
X
(uk − β̃0 )k−i+1
k=ĩ

(k − i + 1)!


aj (k + 1) ,

where
Z

β1

aj (k) = n!
uj̃−1

(β1 − xj̃−1 )j−j̃
(j − j̃)!

Z

xj̃−1

uj̃−2

Z

xk+1

···
uk

dxk ...dxj̃−2 dxj̃−1

j̃−k−1
X ulk+l−1
uj̃−1
β1j−k
= n!
F̄
(
)−
aj (k + 1).
(j − k)! B(j̃−k,j−j̃+1) β1
l!
l=1

A.1.4. Proof of Theorem 3.4. The main idea of the proofs of Theorem
3.4 and 3.5 is as follows. Note that U(i) := D(p(i) ) defined in (15) follow
P
i.i.d.
Γi
the same distribution of Γn+1
, where Γi = ij=1 εj , εj ∼ Exp(1), so that
Γi ∼ Gamma(i, 1). Thus we can approximate
k
P {U(k) > D(g( , b)), for all k0 ≤ k ≤ k1 }.
n
k
≈P {Γ(k) > (n + 1)D(g( , b)), for all k0 ≤ k ≤ k1 }.
n
We take advantage of the joint density of (Γk0 , ..., Γk1 ), which is given by
Lemma 2, while Γn+1 can be approximated by n + 1 when n is reasonably
large.
Accordingly, we apply similar calculation as the proof of Theorem 3.1
except using 5 and Lemma 6 instead.

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 27

A.1.5. Proof of Theorem 3.5. Note that
P (Γk > a + λk, 1 ≤ k ≤ k1 )
Z
Z zk
Z +∞
1
...
e−zk1
=
a+λk1

a+λ(k1 −1)

!

z2

dz1 ...dzk1 −1

dzk1

a+λ

By lemma 4,
Z +∞
(zk1 − a − (k1 − 1)λ) (zk1 − a)k1 −2
=
e−zk1
dzk1
(k1 − 1)!
a+λk1

Z +∞ 
k1 −1
k1 −2
−zk1 (zk1 − a)
−zk1 (k1 − 1)λ(zk1 − a)
=
e
dzk1
−e
(k1 − 1)!
(k1 − 1)!
a+λk1
Let x = zk1 − a,
Z +∞

Z +∞
k1 −1
k1 −2
−a
−x x
−x x
=e
e
dx − λ
e
dx
(k1 − 1)!
(k1 − 2)!
λk1
λk1
=e−a [1 − P (Γk1 ≤ λk1 ) − λ (1 − P (Γk1 −1 ≤ λk1 ))]
=e−a (1 − λ + λFΓ(k1 −1) (λk1 ) − FΓ(k1 ) (λk1 ))
=e−a (1 − λ + hk1 (λ)).
Thus Theorem 3.5 is proved by combining this equation and Lemma 5.
The idea of the proof of Theorem 3.6 is motivated by [14]. Instead of
directly considering the distribution function, we look at the right-tail probability which can be decomposed into the union of disjoint sets.
A.1.6. Proof of Theorem 3.6 and Corollary 3. Let event An,k be defined
S 1
as in Lemma 7. They are disjoint and {Sn ≥ b} = kk=k
An,k . In this proof
0
we mainly focus on approximating P (An,k )
Let dk = (n + 1)D(g( nk , b)), d0k = (n + 1)

k
dD(g( n
,b))
,
dx

d0k,max = (n + 1)

Notice that D(g(x, b)) is convex in x, so dk+j > dk +
and Lemma 7, we have for 1 ≤ k ≤ k1 − 1,

0
dk

n

j. From Lemma 6

dkk
(Γj > dk+j , 1 ≤ j ≤ k1 − k)
k!
0
dkk
dk
≤ P (Γj > dk + j, 1 ≤ j ≤ k1 − k)
k!
n
0
0
k
dk −dk
dk
dk
= e (1 −
+ hk1 −k ( ))
k!
n
n
0
0
dk
dk
= (1 −
+ hk1 −k ( ))fP (dk ) (k).
n
n

P (An,k ) =

dk1 −dk
k1 −k .

28

H. ZHANG ET AL.

Next we need to find the lower bound for P (An,k ).
√
√
We first consider k ≥ k1 − n, here n is chosen for covenience and the
proof works for any nγ , 0 < γ < 1. Note that, due to the convexity, dk+j ≤
0
dk+√n
√
dk + n j, 1 ≤ j ≤ n
0

dk+√n
√
dkk
P (Γj ≥ dk +
j, 1 ≤ j ≤ n)
k!
n
0 √
k
d
d0k+√n
dk −dk
k+ n
√
≥ e (1 −
+ h n(
)
k!
n
n
0
0
d
d
= (1 + o(1))(1 − k + h√n ( k ))fP (dk ) (k).
n
n

P (An,k ) ≥

The last equation is due to Theorem 3.5 and the continuity of D(g(x, b)). We
0
0
√
d
d
can see P (An,k ) → (1 − nk + h√n ( nk ))fP (dk ) (k) uniformly in k, k ≥ k1 − n.
√
We then consider k ≤ k1 − n. In this case, the proof is slightly more
complicated than the first case, however the idea is similar.
Let pn,k (y) = P (Γj > dk+j − dk + y, 1 ≤ j ≤ k1 − k), 0 ≤ y ≤ dk , and fΓk
Rd
be the density function of Γk , then P (An,k ) = 0 k pn,k (y)fΓk (dk − y)dy.
Similar to the proof in the first case,
pn,k (y) ≥ P (Γj ≥ y +

−y

≥e

(1 −

d0k+√n
n

d0k+√n
n

+

= (1 + o(1))e−y (1 −

,1 ≤ j ≤

√

n) − P (

k[
1 −k
√

{Γj ≤ y + dk+j − dk })

j= n

h√

n(

d0k+√n
n

)) −

kX
1 −k
√
j= n

0
Γj
y dk,max
P(
≤ +
)
j
j
n

kX
0
1 −k
d0k
d0
Γj
y dk,max
P(
+ h√n ( k )) −
≤ +
).
n
n
j
j
n
√
j= n

To prove the residual uniformly in k converges to 0, we need to apply
Lemma 3. For y of some lower order, say O(log n)
0
0
y dk,max
log n dk,max
+
≤ √ +
j
n
n
n
log n n + 1
dD(g(x, b))
< √ +
sup {
}
n k0 ≤x≤ k1
dx
n
n

≤δ<1

n

when n is large.

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 29

Then, by Lemma 3, we have the residual
resk ≤

kX
1 −k

√

e

−jI(δ)

√

j= n

e−I(δ) n
→ 0 exponentially uniformly in k.
≤
1 − e−I(δ)

From
(1 + o(1))e−y (1 −

d0k
d0
d0
d0
+ h√n ( k )) − resk ≤ pn,k ≤ e−y (1 − k + h√n ( k )),
n
n
n
n

we can conclude that
d0
d0k
+ h√n ( k )) uniformly in k,
n
n √
n
−I(δ)
e
with error bound resk ≤
.
1 − e−I(δ)

pn,k (y) → e−y (1 −

Let α denote log n for simplicity, α > 1, then, for some constant c > 1,
(16)
Z min{dk ,cα}
d0
d0
e−y (1 − k + h√n ( k ))fΓk (dk − y)dy
Iˆn,k =
n
n
0
√
Z min{dk ,cα}
ne−I(δ) n
→In,k =
pn,k (y)fΓk (dk − y)dy uniformly with error bound ≤
.
1 − e−I(δ)
0
When dk > cα, dk − y < dk < k − 1, fΓk (dk − y) is decreasing in y
P (An,k ) − In,k
(dk − cα)fΓk (dk − cα)
≤
In,k
αfΓk (dk − α)
(dk − cα)ecα−dk (dk − cα)k−1
αeα−dk (dk − α)k−1
(dk − cα) (c−1)α
(c − 1)α k−1
=
e
(1 −
)
α
dk − α
(dk − cα) (c−1)α (k−1) log(1− (c−1)α
)
dk −α
=
e
e
α
(c−1)α
dk
−(k−1)( d −α )
k
≤ eα(c−1) e
α
dk α(c−1)(1− k−1
)
dk
≤ e
α
k
Let x = , y = D(g(x, b)), then
n
nx−1
(n + 1)y −(c−1)( (n+1)y
−1)
=
n
.
α
=

30

H. ZHANG ET AL.

Notice that

nx−1
(n+1)y

is c = sup 2 ≤x≤ k1
n

n

nx−1
> 1, we need (c − 1)( (n+1)y
− 1) > 1 for all x, that

nx−1
nx−1−(n+1)D(g(x,b))

> 1. Therefore

P (An,k )−In,k
In,k

≤

y −γ
αn ,

γ > 0. For k = 1, the result follows d1 → 0,
(17)

In,k → P (An,k ) uniformly in k as n → ∞.

Define P̂ (An,k ) =

R dk
0

e−y (1 −

P̂ (An,k ) − Iˆn,k
P̂ (An,k )

d0k
n

d0

+ h√n ( nk ))fΓk (dk − y)dy,
(dk −y)k−1
(k−1)! dy
R dk (dk −y)k−1
0
(k−1)! dy
dk − (dk − cα)k
)
( k
dkk
k

R min{dk ,cα}

=1−
=1−

0


cα
= 1−
dk


k
≤ exp −c log n .
dk
Since c, dkk > 1, we have
(18)

Iˆn,k → P̂ (An,k ) uniformly in k as n → ∞.
k−1

−y)
Finally, notice that e−y fΓk (dk − y) = e−dk (dk(k−1)!
, then

(19)

d0
d0
P̂ (An,k ) = (1 − k + h√n ( k ))e−dk
n
n
0

Z

dk

0

(dk − y)k−1
dy
(k − 1)!

0

d
d
= (1 − k + h√n ( k ))fP (dk ) (k).
n
n
Combine lemma 5 and equations (16), (17), (18), and (19), we have
P (Sn ≥ b) =(1 + o(1))P (Γk ≤ dk , for some k0 ≤ k ≤ k1 )
=(1 + o(1))

k1
X

P̂ (An,k )

k=k0

=(1 + o(1))

k1
X
k=k0

√
where k ∗ = min{k1 − k, n}.

(1 −

d0k
d0
+ hk∗ ( k ))fP (dk ) (k),
n
n

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 31

To proof Corollary 3, by Theorem 3.6, we have

k1 
X
n+1 0 k
n+1 0 k
P (HC 2004 ≥ b) = (1+o(1))
1−
g ( , b0 ) + hk∗ (
g ( , b0 )) fP ((n+1)g( k ,b0 )) (k) .
n
n
n
n
n
k=k0

A.2. Fundmental Lemmas and Proofs.
Lemma 1. Let U(i) be the ith order statistic of i.i.d. samples from Uniform(0,1), 1 ≤ k0 < k1 ≤ n. Then the joint density of (U(k0 ) , ..., U(k1 ) ) is
f(U(k

0)

,...,U(k1 ) ) (zk0 , ..., zk1 )

=

n!
z k0 −1 (1 − zk1 )n−k1 , 0 < zk0 < ... < zk1 < 1.
(n − k1 )!(k0 − 1)! k0

Proof. Standard result from formula (2.2.2) of the book Order Statistics[36].
Lemma
Pk2. Let εi be i.i.d. exponential random variables with parameter
1. Γk = i=1 εi , 1 ≤ k0 < k1 ≤ n. Then the joint density of (Γk0 , ..., Γk1 )
f (zk0 , ..., zk1 ) =

zkk00 −1
(k0 − 1)!

e−zk1 , 0 < zk0 < ... < zk1 < ∞.

Specially, when k0 = 1, the joint density is
f (z1 , ..., zk1 ) = e−zk1 , 0 < z1 < ... < zk1 < ∞.
Proof. The proof follows the deduction process shown in Theorem 1.1
by Mathai and Moschopoulos [37]. Note that their theorem only considers
the case of k0 = 1, but the deduction idea can be applied to general k0 .
P
Lemma 3. Let c < 1, I(c) = c − 1 − log(c) > 0, Γn = ni=1 i , where i
are i.i.d exponential distributed with parameter 1, then
P(

Γn
≤ c) ≤ e−nI(c) .
n

Proof. Let i be exp(1) distributed. The log moment generating function
of −i is Λ(t) = − log(1 + t), t > −1. Then the convex rate function is
Λ∗ (x) = sup tx − Λ(t)
t>−1

= −x − 1 − log(−x).
By Cramer’s theorem, let x = −c, we have desired result.

32

H. ZHANG ET AL.

Lemma 4.
Z

(Abel-Goncharov Polynomial)

zk 1

Z

z2

dz1 ...dzk1 −1 =

...

a+λ(k1 −1)

a+λ

(zk1 − a − λ(k1 − 1)) (zk1 − a)k1 −2
.
(k1 − 1)!

Proof. We will prove by induction.
When k1 = 2, it’s easily shown that
Z z2
(z2 − a − λ(2 − 1)) (z2 − a)2−2
dz1 = z2 − a − λ =
.
(2 − 1)!
a+λ
Assume it’s true for k1 = k,
Z zk
Z z2
(zk − a − λ(k − 1)) (zk − a)k−2
.
...
dz1 ...dzk−1 =
(k − 1)!
a+λ(k−1)
a+λ
Then for k1 = k + 1,
Z zk+1 Z

zk

Z

z2

...
a+λk
Z zk+1

=
a+λk
zk+1

Z
=

a+λk

a+λ(k−1)

dz1 ...dzk−1 dzk
a+λ

(zk − a − λ(k − 1)) (zk − a)k−2
dzk
(k − 1)!
Z zk+1
(zk − a)k−1
λ(zk − a)k−2
dk −
dk
(k − 1)!
(k − 2)!
a+λk

(λk)k
λ(zk+1 − a)k−1 λ(λk)k−1
(zk+1 − a)k
−
−
+
=
k!
k!
(k − 1)!
(k − 1)!
k−1
(zk+1 − a) (zk+1 − a − λk)
.
=
k!
This finished the proof.
Lemma 5. Let c = (ck0 , ..., ck1 ), 1 ≤ k0 < k1 ≤ n be a finite sequence of
increasing numbers. Then given  > 0, as n → ∞
P(

Γj
> cj , k0 ≤ j ≤ k1 ) = (1 + O())P (Γj > (n + 1)cj , k0 ≤ j ≤ k1 ).
Γn+1

Proof. Given  > 0, let L(c; n, ) and R(c; n, ) be the lower bound and

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 33
Γ

j
> cj , k0 ≤ j ≤ k1 ),
upper bound for P ( Γn+1

L(c; n, ) =P (Γj > (1 + )(n + 1)cj , k0 ≤ j ≤ k1 )
Γn+1
−P (
∈ [1 − , 1 + ]c ),
n+1
R(c; n, ) =P (Γj > (1 − )(n + 1)cj , k0 ≤ j ≤ k1 )
Γn+1
∈ [1 − , 1 + ]c ).
+P (
n+1
By Lemma 3,
P(

Γn+1
1
∈ [1 − , 1 + ]c ) ≤ 2 e−nΛ(1−) → 0 exponentially.
n+1


Because of the continuity of the multivariate gamma distribtution,
L(c; n, ) →(1 + O())P (Γj > (n + 1)cj , k0 ≤ j ≤ k1 ),
R(c; n, ) →(1 + O())P (Γj > (n + 1)cj , k0 ≤ j ≤ k1 ).

Lemma 6.
Let (d1 , ..., dk1 ) be a sequence of nondecreasing and nonnegative numbers. Qj (dj+1 , ..., dk1 ) = P {Γk ≥ dk+j , 1 ≤ k ≤ k1 − j},
0 ≤ j ≤ k1 − 1. F̄Γk (x) is the survival function of Gamma distribution
with shape parameter k and scale parameter 1. Then for l = 2, 3, ..., k1
Qk1 −l = F̄Γl (dk1 ) −

l−1 j
X
dk

1 −l+j

j=1

j!

Qk1 −l+j .

with Qk1 −1 = F̄Γ1 (dk1 ), and, for k0 ≥ 1, the joint survival probability
P {Γk ≥ dk , k0 ≤ k ≤ k1 } = F̄Γk1 (dk1 ) −

kX
1 −1
j=k0

djj
j!

Qj .

34

H. ZHANG ET AL.

Proof.
Q0 = P {Γk ≥ dk , 1 ≤ k ≤ k1 }
Z z2
Z zk
Z +∞
1
−zk1
dz1 ...dzk1 −1 dzk1
...
e
=
+∞

=

e−zk1

Z

+∞

e

=

−zk1

Z

+∞

=

e

−zk1

Z

+∞

=
dk1

e−zk1

zk 1

Z

z2 dz2 ...dzk1 −1 dzk1 − d1 Q1
d2

Z
...

djj
j!

zk0 +1

dk0

(k1 − 1)!

j=1

z3

...

zk 1

kX
1 −1

(z2 − d1 )dz2 ...dzk1 −1 dzk1
d2

zkk11 −1

= F̄Γk1 (dk1 ) −

z3

...

dk1 −1

dk1

Z

Z

dk1 −1

dk1

Z

zk 1

dk1 −1

dk1

Z

d1

dk1 −1

dk1

Z

dzk1 −

zkk00 −1
(k0 − 1)!
kX
1 −1
j=1

djj
j!

dzk0 ...dzk1 −1 dzk1 −

kX
0 −1
j=1

djj
j!

Qj

Qj

Qj .

The rest of recursive formulae can be derived similarly and from the last
third equation.
Also notice that {Γk ≥ dk , k0 ≤ k ≤ k1 } = {Γk ≥ dk , 1 ≤
S
k ≤ k1 } {some Γk < dk , 1 ≤ k ≤ k0 − 1 and Γk ≥ dk , k0 ≤ k ≤ k1 }, by
Lemma 7 we have
P {Γk ≥ dk , k0 ≤ k ≤ k1 } = Q0 +

kX
0 −1
j=1

djj
j!

= F̄Γk1 (dk1 ) −

Qj

kX
1 −1
j=k0

djj
j!

Qj .

Lemma 7. Let (d1 , ..., dk1 ) be a sequence of increasing numbers, An,k =
{Γk ≤ dk , Γk+l > dk+l , 1 ≤ l ≤ k1 − k}, 1 ≤ k ≤ k1 − 1. Then
P (An,k ) =

dkk
Qk .
k!

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 35

Proof.
P (An,k ) = P {Γk ≤ dk , Γk+l > dk+l for all 1 ≤ l ≤ k1 − k}
Z zk+1 Z dk
Z zk
Z +∞
1
zkk−1
...
e−zk1
=
dzk dzk+1 ...dzk1 −1 dzk1
(k − 1)!
0
dk+1
dk1 −1
dk1
Z zk+1 k
Z zk
Z +∞
1
dk
−zk1
...
e
=
dzk+1 ...dzk1 −1 dzk1
dk+1 k!
dk −1
dk
1

1

dk
= k P {Γl ≥ dk+l , 1 ≤ l ≤ k1 − k}
k!
dk
= k Qk .
k!

REFERENCES
[1] D. Donoho and J. Jin, “Higher criticism for detecting sparse heterogeneous mixtures,”
The Annals of Statistics, vol. 32, no. 3, pp. 962–994, 2004.
[2] D. Donoho and J. Jin, “Higher criticism thresholding: Optimal feature selection when
useful features are rare and weak,” Proceedings of the National Academy of Sciences
of the United States of America, vol. 105, pp. 14790–14795, Sep 30 2008.
[3] D. Donoho and J. Jin, “Higher criticism for large-scale inference: especially for rare
and weak effects,” Statistical Science, vol. 30, no. 1, pp. 1–25 DOI: 10.1214/14–
STS506, 2015.
[4] L. Jager and J. A. Wellner, “Goodness-of-fit tests via phi-divergences,” The Annals
of Statistics, pp. 2018–2053, 2007.
[5] M. Noé, “The calculation of distributions of two-sided kolmogorov-smirnov type
statistics,” The Annals of Mathematical Statistics, pp. 58–64, 1972.
[6] M. Noé and G. Vandewiele, “The calculation of distributions of kolmogorov-smirnov
type statistics including a table of significance points for a particular case,” The
Annals of Mathematical Statistics, vol. 39, no. 1, pp. 233–241, 1968.
[7] V. Kotel’Nikova and E. Chmaladze, “On computing the probability of an empirical
process not crossing a curvilinear boundary,” Theory of Probability & Its Applications,
vol. 27, no. 3, pp. 640–648, 1983.
[8] G. R. Shorack and J. A. Wellner, Empirical processes with applications to statistics,
vol. 59. SIAM, 2009.
[9] G. Steck, “The smirnov two sample tests as rank tests,” The Annals of Mathematical
Statistics, pp. 1449–1466, 1969.
[10] M. Breth, “On a recurrence of steck,” Journal of Applied Probability, pp. 823–825,
1976.
[11] H. Ruben, “On the evaluation of steck’s determinant for rectangle probabilities of
uniform order statistics,” Communications in Statistics-Theory and Methods, vol. 5,
no. 6, pp. 535–543, 1976.
[12] I. J. Barnett and X. Lin, “Analytical p-value calculation for the higher criticism test
in finite-d problems,” Biometrika, vol. 101, no. 4, pp. 964–970, 2014.

36

H. ZHANG ET AL.

[13] A. Moscovich, B. Nadler, C. Spiegelman, et al., “On the exact berk-jones statistics and
their p-value calculation,” Electronic Journal of Statistics, vol. 10, no. 2, pp. 2329–
2354, 2016.
[14] J. Li and D. Siegmund, “Higher criticism: p-values and criticism,” Annals of Statistics,
vol. 43, no. 3, pp. 1323–1350, 2015.
[15] F. Eicker, “The asymptotic distribution of the suprema of the standardized empirical
processes,” The Annals of Statistics, pp. 116–138, 1979.
[16] D. Jaeschke, “The asymptotic distribution of the supremum of the standardized empirical distribution function on subintervals,” The Annals of Statistics, vol. 7, no. 1,
pp. 108–115, 1979.
[17] Y. I. Ingster, “Some problems of hypothesis testing leading to infinitely divisible
distributions,” Mathematical Methods of Statistics, vol. 6, no. 1, pp. 47–69, 1997.
[18] Y. I. Ingster, “Minimax detection of a signal for in -balls,” Mathematical Methods of
Statistics, vol. 7, no. 4, pp. 401–428, 1998.
[19] J. Tukey, “The higher criticism.” Course Notes, Statistics 411, Princeton University.,
1976.
[20] E. Arias-Castro, E. J. Candès, and Y. Plan, “Global testing under sparse alternatives:
Anova, multiple comparisons and the higher criticism,” The Annals of Statistics,
vol. 39, no. 5, pp. 2533–2556, 2011.
[21] J. Shao, Mathematical Statistics. Springer Verlag, 2010.
[22] R. H. Berk and D. H. Jones, “Goodness-of-fit test statistics that dominate the kolmogorov statistics,” Probability Theory and Related Fields, vol. 47, no. 1, pp. 47–59,
1979.
[23] D. B. Goldstein, “Common genetic variation and human traits,” New England Journal of Medicine, vol. 360, no. 17, pp. 1696–1698, 2009.
[24] A. L. Price, N. J. Patterson, R. M. Plenge, M. E. Weinblatt, N. A. Shadick, and
D. Reich, “Principal components analysis corrects for stratification in genome-wide
association studies,” Nature genetics, vol. 38, no. 8, pp. 904–909, 2006.
[25] P. McCullagh and J. A. Nelder, Generalized linear models, vol. 37. CRC press, 1989.
[26] D. J. Schaid, C. M. Rowland, D. E. Tines, R. M. Jacobson, and G. A. Poland, “Score
tests for association between traits and haplotypes when linkage phase is ambiguous,”
The American Journal of Human Genetics, vol. 70, no. 2, pp. 425–434, 2002.
[27] R. Duerr, K. Taylor, S. Brant, J. Rioux, M. Silverberg, M. Daly, A. Steinhart,
C. Abraham, M. Regueiro, A. Griffiths, et al., “A genome–wide association study
identifies il23r as an inflammatory bowel disease gene,” Science Signalling, vol. 314,
no. 5804, p. 1461, 2006.
[28] J.-P. Hugot, M. Chamaillard, H. Zouali, S. Lesage, J.-P. Cézard, J. Belaiche, S. Almer,
C. Tysk, C. A. O’Morain, M. Gassull, et al., “Association of nod2 leucine-rich repeat
variants with susceptibility to crohn’s disease,” Nature, vol. 411, no. 6837, pp. 599–
603, 2001.
[29] Y. Ogura, D. K. Bonen, N. Inohara, D. L. Nicolae, F. F. Chen, R. Ramos, H. Britton,
T. Moran, R. Karaliuskas, R. H. Duerr, et al., “A frameshift mutation in nod2 associated with susceptibility to crohn’s disease,” Nature, vol. 411, no. 6837, pp. 603–606,
2001.
[30] D. Maglott, J. Ostell, K. D. Pruitt, and T. Tatusova, “Entrez gene: gene-centered
information at ncbi,” Nucleic acids research, vol. 39, no. suppl 1, pp. D52–D57, 2011.
[31] P. Chamouard, Z. Richert, N. Meyer, G. Rahmi, and R. Baumann, “Diagnostic value
of c-reactive protein for predicting activity level of crohns disease,” Clinical Gastroenterology and Hepatology, vol. 4, no. 7, pp. 882–887, 2006.
[32] G. Trikudanathan, P. G. Venkatesh, and U. Navaneethan, “Diagnosis and therapeu-

DISTRIBUTIONS OF SIGNAL-DETECTION METHODS IN FINITE CASES 37

[33]
[34]

[35]

[36]
[37]

tic management of extra-intestinal manifestations of inflammatory bowel disease,”
Drugs, vol. 72, no. 18, pp. 2333–2349, 2012.
W. D. James, T. Berger, and D. Elston, Andrew’s diseases of the skin: clinical dermatology. Elsevier Health Sciences, 2011.
S. Yuvaraj, S. H. Al-Lahham, R. Somasundaram, P. A. Figaroa, M. P. Peppelenbosch,
and N. A. Bos, “E. coli-produced bmp-2 as a chemopreventive strategy for colon
cancer: a proof-of-concept study,” Gastroenterology research and practice, vol. 2012,
2012.
M. L. Slattery, A. Lundgreen, J. S. Herrick, S. Kadlubar, B. J. Caan, J. D. Potter,
and R. K. Wolff, “Genetic variation in bone morphogenetic protein and colon and
rectal cancer,” International Journal of Cancer, vol. 130, no. 3, pp. 653–664, 2012.
H. David and H. Nagaraja, “Order statistics,” 2003.
A. M. Mathai and P. G. Moschopoulos, “On a multivariate gamma,” Journal of
Multivariate Analysis, vol. 39, no. 1, pp. 135–153, 1991.
Department of Mathematical Sciences and
Program of Bioinformatics and Computational Biology and
Program of Data Science
100 Institute Road
Worcester, MA 01609
E-mail: zheyangwu@wpi.edu
URL: https://www.wpi.edu/people/faculty/zheyangwu

