arXiv: arXiv:0000.0000

Permutation p-value approximation via
generalized Stolarsky invariance

arXiv:1603.02757v2 [] 9 Aug 2017

Hera Y. He, Kinjal Basu, Qingyuan Zhao, Art B. Owen
Department of Statistics
Sequoia Hall
Stanford University
Stanford, CA 94305
e-mail: (yhe1, kinjal, qyzhao, owen)@stanford.edu
Abstract: It is common for genomic data analysis to use p-values from a
large number of permutation tests. The multiplicity of tests may require
very tiny p-values in order to reject any null hypotheses and the common practice of using randomly sampled permutations then becomes very
expensive. We propose an inexpensive approximation to p-values for two
sample linear test statistics, derived from Stolarsky’s invariance principle.
The method creates a geometrically derived set of approximate p-values
for each hypothesis. The average of that set is used as a point estimate p̂
and our generalization of the invariance principle allows us to compute the
variance of the p-values in that set. We find that in cases where the point
estimate is small the variance is a modest multiple of the square of the point
estimate, yielding a relative error property similar to that of saddlepoint
approximations. On a Parkinson’s disease data set, the new approximation
is faster and more accurate than the saddlepoint approximation. We also
obtain a simple probabilistic explanation of Stolarsky’s invariance principle.

1. Introduction
Permutation methods are commonly used to obtain p-values in genomic applications. They make only modest assumptions and they have a direct intuitive
interpretation that appeals to biologists in collaborations. In even modestly
large data sets, the exact permutation p-value becomes too expensive to compute. Then Monte Carlo sampling of random permutations becomes a standard
approach. Genomic applications commonly require thousands or more of hypotheses to be tested, and then multiplicity adjustment requires that some small
p-values be obtained if any null hypotheses are to be rejected. When p-values
below  are required to reject H0 , then Knijnenburg et al. (2009) recommend
doing at least 10/ random permutations. As a result, even Monte Carlo sampling for permutation tests can be prohibitively expensive, and hence it pays to
search for fast approximations to the permutation p-value.
In this paper we develop rapidly computable approximations to some permutation p-values. The p-values we consider are for a difference in group means.
The approximations are based on ideas from spherical geometry and discrepancy,
related to the Stolarsky invariance principle (Stolarsky, 1973). As described below, the resulting approximations prove to be very accurate for the tiny p-values
where permutation methods are most difficult to use.
1

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

2

We begin with some background on the genomic motivation of our work.
Then we transition to spherical geometry.
Genomic context
The specific problem that motivated us is testing for sets of genes associated
with Parkinson’s disease (Larson and Owen, 2015). More details about this work
are given in the the first author’s dissertation (He, 2016).
In these data sets, there are m0 subjects without Parkinson’s disease and m1
subjects with it. One can test whether Parkinson’s disease is associated with
an individual gene by doing a t-test comparing gene expression levels in tissue
samples from the two groups of subjects. Biological interest is often summarized
more by gene sets rather than individual genes. Gene sets have two advantages:
small but consistent associations of many genes with the test condition can
raise power, and, the gene sets themselves often connect better to biological
understanding than do individual genes. The null hypothesis is that the subject
condition does not affect expression levels of any gene in the gene set.
There is a very large literature on testing for significant associations between
a condition and the genes in a gene set. See Ackermann and Strimmer (2009)
for an overview of the main concepts and methods in that literature. They did
an extensive comparison of 261 different gene set testing methods and identified
two families of winning methods. Let tg be the ordinary two sample t statistic
for comparing the average expression level of gene g between two conditions and
let G be a set ofPgenes of interest. They
P found that linear and quadratic test
statistics LG = g∈G tg and QG = g∈G t2g had the best power, along with
some simple approximations to those two statistics. These methods performed
better than some subtantially more complicated proposals. The statistic LG was
proposed by Jiang and Gentleman (2007). The statistic LG and similar ones did
best when expression differences between the two conditions tend to have the
same sign, for each g ∈ G. If large but oppositely signed treatment effects occur,
then QG and approximations to it do best.
The genes in a gene set are ordinarily correlated with each other, even if
they are independent of the treatment condition. The correlations makes it
difficult to find the null distributions of LG and QG , even with parametric model
assumptions. In a permutation
 analysis, like Ackermann and Strimmer (2009)
use, we consider all N = mn1 different ways to select a subset π containing m1
of the n = m0 + m1 subjects. Let QπG be the test statistic recomputed as if those
m1 subjects had been the affected group. Then the permutation p-value for QG
is
1 X
1QπG >QG .
p = pG =
N π
Note that the smallest possible value for p is 1/N .
When N is too large for a permutation test to be computationally feasible,
a standard practice is to estimate p via randomly sampled permutations of the
treatment label as proposed by Barnard (1963). For ` = 1, . . . , M − 1 we let

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

3

π(`) be the affected group after a randomization of the treatment labels. We let
π(0) be the original allocation. Then the Monte Carlo estimate
p̂ =

M −1
1 X
1Qπ(`) >QG
G
M
`=0

is used as an estimate of p, for the quadratic statistic. In this Monte Carlo,
the true permutation p-value p is the unknown parameter and p̂ is the sample
estimate of p. Note that p̂ > 1/M because we have included the original allocation in the numerator. Failure to include the original allocation π(0) can lead
to p̂ = 0 which is very undesirable. We call 1/M the granularity limit. When
p is quite small, an enormous number M of simulations may be required to
get an accurate estimate of it. For instance, in genome wide association studies (GWAS) the customary threshold for significance is  = 5 × 10−8 , making
permutation methods prohibitively expensive, or even infeasible. For a recent
discussion of p-value thresholds in GWAS, see Fadista et al. (2016).
In this paper, we work with one of Ackermann and Strimmer’s (2009) approximations to LG . Let Xi = 1 if subject i is in condition 1 and Xi = 0 for
condition 0. Let Yig be the expression level of gene g for subjectq
i. Let ρ̂g be
√
the sample correlation between Xi and Yig . Then tg = n − 2ρ̂g / 1 − ρ̂2g and
. √
a first order Taylor approximation gives tg = n − 2(ρ̂g + ρ̂3g /2). When many
small correlations ρ̂g contribute to the signal, then summing ρ̂g gives a test
statistic that is almost equivalent to summing tg . Ackermann and Strimmer
(2009) found that
n
X 1X
Xi − X̄ Ygi − Ȳg
n i=1 sX
sg

(1.1)

g∈G

was in the same winning set of methods as LG , where sX and sg are standard deviations of Xi and Ygi respectively. They also considered using pooled
variance estimates in place of sg but found no advantage to doing so, perhaps
because
min(m0 , m1 ) was at least 10 in their simulations. Letting Yi = YGi ≡
P
g∈G Ygi /sg , we may rewrite (1.1) as
n
X
i=1

Yi − Ȳ
Xi − X̄
kX − 1n X̄k kY − 1n Ȳ k

(1.2)

multiplied by a constant that only depends on n and hence does not affect p.
Equation (1.2) describes a test statistic that is a plain Euclidean inner product
of two unit vectors in Rn . Here x0 has i’th component (Xi − X̄)/kX − 1n X̄k
and y 0 is similar.

There are N = mn1 distinct vectors found by permuting the entries in x. We
label them x0 , x1 , . . . , xN −1 with x0 being the original one. Letting ρ̂ = xT
0 y0
we find that one and two-sided p-values for a linear statistic are
N −1
1 X
1xTk y0 >ρ̂
N
k=0

and

N −1
1 X
1|xTk y0 |>|ρ̂|
N
k=0

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

4

repectively. We prefer two-sided test statistics, but we will study one-sided ones
first and then translate our results to two-sided ones.
Spherical geometry
We are now ready to make a geometric interpretation. Let Sd = {z ∈ Rd+1 |
z T z = 1} be the d-dimensional unit sphere. Our data x0 , y 0 are in a subset
of Sn−1 orthogonal to 1n . That subset is isomorphic to Sn−2 and so we work
mostly with d = n − 2.
Given a point y ∈ Sd , the points z that are closest to y comprise a spherical
cap. The spherical cap of center y and height t ∈ [−1, 1] is C(y; t) = {z ∈ Sd |
hy, zi > t}. By symmetry, z ∈ C(y; t) if and only if y ∈ C(z; t). The one-sided
linear p-value is the fraction of xk for 0 6 k < N that belong to C(y 0 ; ρ̂). A
natural, but crude approximation to p is then
p̂1 (ρ̂),

where

p̂1 (t) ≡

vol(C(y 0 ; t))
.
vol(Sd )

Stolarsky’s invariance principal gives a remarkable description of the accuracy
of this approximation p̂1 . The squared L2 spherical cap discrepancy of points
x0 , x1 , . . . , xN −1 ∈ Sd is
Z 1Z
2
L2 (x0 , . . . , xN −1 ) =
|p̂1 (t) − p(z, t)|2 dσd (z) dt
−1

Sd

PN −1
where p(z, t) = (1/N ) k=0 1xk ∈C(z,t) and σd is the uniform (Haar) measure
on Sd . Stolarsky (1973) shows that
dωd
× L2 (·)2 =
ωd+1

Z

Z
kx − yk dσd (x) dσd (y) −

Sd

Sd

N −1
1 X
kxk − xl k
N2

(1.3)

k,l=0

where ωd is the (surface) volume of Sd . Equation (1.3) relates the mean squared
error of p̂1 to the mean absolute Euclidean distance among the N points. In our
applications, the N points will be the distinct permuted values of x0 , but (1.3)
holds for an arbitrary set of N points xk .
The left side of (1.3) is, up to normalization, a mean squared discrepancy over
spherical caps. This average of (p̂1 − p)2 includes p-values of all sizes between 0
and 1. It is not then a very good accuracy measure when p̂1 (ρ̂) turns out to be
very small, such as 10−6 . It would be more useful to get such a mean squared
error taken over caps of exactly the size p̂1 (ρ̂), and no others.
Brauchart and Dick (2013) consider quasi-Monte Carlo (QMC) sampling in
the sphere. They generalize Stolarsky’s discrepancy formula to include a weighting function on the height t. By specializing their formula, we get an expression
for the mean of (p̂1 − p)2 over spherical caps of any fixed size.
Discrepancy theory plays a prominent role in QMC (Niederreiter, 1992),
which is about approximating an integral by a sample average. The present

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

5

setting is a reversal of QMC: the discrete average p over permutations is the
exact value we seek, and the integral over a continuum is the approximation p̂.
A second difference is that the QMC literature focusses on choosing N points
to minimize a criterion such as (1.3), whereas here the N points are determined
by the problem.
As we will show below, the estimate p̂1 is the average of p over all spherical
caps C(y; ρ̂) under a uniform distribution, i.e., y ∼ U(Sd ). Those caps have the
same volume as C(y 0 ; ρ̂).
In addition to specializing from caps C(y; t) with t = ρ̂ we can also specialize
to caps whose centers y more closely resemble y 0 . Suppose that
y 0 ∈ Y ⊂ Sd .
Then we know that p(y 0 , ρ̂) 6 supy∈Y p(y, ρ̂), which is a conservative permutation p-value. We are generally unable to compute this quantity but in some
instances we can form a reference distribution y ∼ U(Y) and compute both
E(p(y, ρ̂) | y ∈ Y) and Var(p(y, ρ̂) | y ∈ Y), the mean and variance of p(y, ρ̂)
under this distribution.
The simplest reference distribution we use has Y1 = Sd . We have found
that the set Y2 = {y ∈ Sd | y T x0 = y T
0 x0 } and some generalizations yield
especially useful reference distributions. Generalizations of the Stolarsky formula
allow us to compute Var(p(y, ρ̂) | y ∈ Y2 ). In some of our numerical results
from Section 7, we find that Var(p(y, ρ̂) | y ∈ Y2 ) is so small that the true
permutation p-value p(y 0 , ρ̂) must be of the same order of magnitude as the
estimate p̂2 ≡ E(p(y, ρ̂) | y ∈ Y2 ) that we study at length.
We obtain p̂2 and the mean square discrepancy over its reference distribution by further extending Brauchart and Dick’s generalization of Stolarsky’s
invariance. More generally, we can replace the constraint hy, x0 i = hy 0 , x0 i by
hy, xc i = hy 0 , xc i for any individual c ∈ {0, 1, . . . , N − 1}. Our estimate p̂3 takes
xc to be whichever permuted point xk happens to be closest to y 0 .
Smaller sets Y could be even better than Y2 . In the extreme, if we could work
with Y that satisfies hy, xk i = hy 0 , xk i, for all 0 6 k < N then all y ∈ Y would
have p(y, ρ̂) = p(y 0 , ρ̂) = p and the mean over y ∈ Y would have no error. Any
smaller set is only useful if we can efficiently compute with it.
Although we found these results via invariance, we can also obtain them via
probabilistic arguments. As a consequence we have a probabilistic derivation
of Stolarsky’s formula. Bilyk et al. (2016) have independently found this connection. Some of our results are for arbitrary x, but our best computational
formulas are for the case where the variable x is binary, as it is for the Parkinson’s disease data sets.
Outline
The rest of the paper is organized as follows. Section 2 presents some context
on permutation tests and gives some results from spherical geometry. In Section 3 we use Stolarsky’s invariance principle as generalized by Brauchart and

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

6

Dick (2013) to obtain the mean squared error between the true p-value and its
continuous approximation p̂1 , averaging over all spherical caps of volume p̂1 .
This section also has a probabilistic derivation of that mean squared error. In
Section 4 we describe some finer approximations p̃ for the p-value. These use
the set Y2 to condition on not just the volume of the spherical cap but also on
its distance from the original data point x0 , or from some other point, such as
the closest permutation of x0 to y 0 . By always including the original point we
ensure that p̃ > 1/N . That is a desirable property because the true permutation p-value cannot be smaller than 1/N . In Section 5 we modify the proof in
Brauchart and Dick (2013), to further generalize their invariance results to include the mean squared error of the finer approximations. Section 6 extends our
estimates to two-sided testing. Section 7 illustrates our p-value approximations
numerically. We see that an RMS error in the estimate p̂2 is of the same order of
magnitude as p̂2 itself. That is, p̂2 has a relative error property like saddlepoint
estimates do. Section 8 makes a numerical comparison to saddlepoint methods
in simulated data. The saddlepoint estimates come out more accurate than p̂2
but are biased low in the simulated examples. Section 9 compares the accuracy
of our approximations to each other and to the saddlepoint approximation for
6180 gene sets and some Parkinson’s disease data sets. In the data examples,
the new approximations come out closer to some gold standard estimates (based
on large Monte Carlo samples) than the saddlepoint estimates do, which once
again are biased low. From Table 6.3 of He (2016), the saddlepoint computations
take roughly 30 times longer than p̂2 does. Section 10 draws some conclusions
and discusses the challenges in getting a computationally feasible p-value that
accounts for both sampling uncertainty of (x, y) and the uncertainty in p̂ as an
estimate of p. Most of the proofs are in the Appendix, Section 11.
Software
The proposed approximations are implemented in the R package pipeGS on
CRAN. Given a binary input label and a gene expression measurement matrix,
it computes our three p value approximations for the linear gene set statistics.
Those statistics, p̂1 , p̂2 , p̂3 , are mentioned above and then presented in more
detail in Section 4. We provide an implementation of the saddlepoint approximation in the package as well.
2. Background and notation
The raw data contain points (Xi , Yi ) for i = 1, . . . , n, where Yi may be a composite quantity derived from all Ygi for g belonging to a gene set G, such
as Yi = YGi just before (1.2). We center and scale vectors (X1 , X2 , . . . , Xn )
and (Y1 , Y2 , . . . , Yn ) yielding x0 , y 0 ∈ Sd for d = n − 1. Both points belong to
{z ∈ Sn−1 | z T 1n = 0}. We can use an orthogonal matrix to rotate the points
of this set onto Sn−2 × {0}. As a result, we may simply work with x0 , y 0 ∈ Sd
where d = n − 2.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

7

The sample correlation of these variables is ρ̂ = xT
0 y 0 = hx0 , y 0 i. We use
hx0 , y 0 i when we find that geometrical thinking is appropriate and to conform
with Brauchart and Dick (2013). We use xT
0 y 0 to emphasize computational or
algebraic connotations.
Here we develop approximations to the one-sided p-value as that simplifies
notation. Section 6 shows how to obtain the corresponding two-sided p-values.
We assume that ρ̂ > 0 for otherwise p̂1 is going to be too large to be interesting.
For instance with m0 = m1 , ρ̂ 6 0 implies that p̂1 > 1/2.
Our proposals are computationally most attractive in the case where Xi takes
on just two values, such as 0 and 1. Then ρ̂ is a two-sample test statistic for
a difference in means. When there are m0 observations with
p Xi = 0 and m1
m1 /(nm0 ) and m1
with Xi = 1 then x0 contains
m
components
equal
to
−
0
p
components equal to + m0 /(nm1 ). Computational costs are often sensitive to
the smaller sample size, m ≡ min(m0 , m1 ).

1
For this two-sample case there are only N = m0m+m
distinct permuta0
tions of x0 . We have called these x0 , x1 , . . . , xN −1 and the true p value is
PN −1
PN −1
p = (1/N ) k=0 1(xT
k y 0 > ρ̂) = (1/N )
k=0 1(xk ∈ C(y 0 ; ρ̂)).
Now suppose that there are exactly r indices for which xk is positive and x`
is negative. There are then r indices with the reverse pattern too. We say that
xk and x` are at ‘swap distance r’ because r zeros from xk swapped positions
with ones to yield x` . In that case we easily find that
u(r) ≡ hxk , x` i = 1 − r

 1
1 
+
.
m0
m1

(2.1)

We need some geometric properties of the unit sphere and spherical caps. The
surface volume of Sd is ωd = 2π (d+1)/2 /Γ((d + 1)/2). We use σd for the volume
element in Sd normalized so that σd (Sd ) = 1. The spherical cap C(y; t) = {z ∈
Sd | z T y > t} has volume
(

1
I1−t2 d2 , 21 ,
06t61
2

σd (C(y; t)) =
d 1
1
2
1 − 2 I1−t 2 , 2 , −1 6 t < 0
where It (a, b) is the incomplete beta function
1
It (a, b) =
B(a, b)

Z

t

xa−1 (1 − x)b−1 dx

0

R1
with B(a, b) = 0 xa−1 (1 − x)b−1 dx. Obviously, this volume is 0 if t < −1 and
it is 1 if t > 1. This volume is independent of y so we may write σd (C(· , t)) for
the volume.
Our first approximation of the p-value is p̂1 (ρ̂) = σd (C(y; ρ̂)). We remarked
earlier that this approximation equates a discrete fraction to a volume ratio. We
show in Proposition 2 that p̂1 = E(p | hx0 , yi = ρ̂) for y ∼ U(Sd ) as y 0 would
if the original Yi were IID Gaussian. In Theorem 4, we find Var(p̂1 ) under this
assumption.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

8

d
d
We frequently
√ need to project y ∈ S onto a point x ∈ S . In this representation y = tx + 1 − t2 y ∗ where t = y T x ∈ [−1, 1] and y ∗ ∈ {z ∈ Sd | z T x = 0}
which is isomorphic to Sd−1 . The coordinates t and y ∗ are unique. From equation
(A.1) in Brauchart and Dick (2013) we get

dσd (y) =

ωd−1
(1 − t2 )d/2−1 dt dσd−1 (y ∗ ).
ωd

(2.2)

In their case x was (0, 0, . . . , 1).
The intersection of two spherical caps of common height t is
C2 (x, y; t) ≡ C(x; t) ∩ C(y; t).
We will need the volume of this intersection. Lee and Kim (2014) give a general
solution for spherical cap intersections without requiring equal heights. They
enumerate 25 cases, but our case does not correspond to any single case of
theirs and so we obtain the formula we need directly, below. We suspect it must
be known already, but we were unable to find it in the literature.
Lemma 1. Let x, y ∈ Sd and −1 6 t 6 1 and put u = xT y. Let V2 (u; t, d) =
σd (C2 (x, y; t)). If u = 1, then V2 (u; t, d) = σd (C(x; t)). If −1 < u < 1, then
Z
d
ωd−1 1
(1 − s2 ) 2 −1 σd−1 (C(y ∗ ; ρ(s))) ds,
V2 (u; t, d) =
ωd t
p
where ρ(s) = (t − su)/ (1 − s2 )(1 − u2 ). Finally, for u = −1,
(
0,
t>0
V2 (u; t, d) = ωd−1 R |t|
−1
2 d
(1 − s ) 2 ds, else.
ωd
−|t|

(2.3)

(2.4)

Proof. Let z ∼ U(Sd ). Then V2 (u; t, d) = σd (C2 (x, y; t)) = Pr(z ∈ C2 (x, y; t)).
If u = 1 then x = y and√so C2 (x, y; t) = C(x; t).√For u < 1, we project y and
z onto x, via z = sx + 1 − s2 z ∗ and y = ux + 1 − u2 y ∗ . Now
Z
V2 (u; t, d) =
1(hx, zi > t)1(hy, zi > t) dσ(z)
Sd
1

Z
=

d
ωd−1
(1 − s2 ) 2 −1
ωd
p
p
1(su + 1 − s2 1 − u2 hy ∗ , z ∗ i > t) dσd−1 (z ∗ ) ds.

1(s > t)
−1

Z
×

Sd−1

If u > −1 then this reduces to (2.3). For u = −1 we get
V2 (u; t, d) =
which reduces to (2.4).

ωd−1
ωd

Z

1

−1

d

1(s > t)1(−s > t)(1 − s2 ) 2 −1 ds.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

9

When we give probabilistic arguments and interpretations we do so for a random center y of a spherical cap. That random center is taken from two reference
distributions. Those are distributions 1 and 2 below. Reference distribution 1
is illustrated in Figure 1. Distribution 2 is illustrated in Figure 2 of Section 4
where we first use it.
Reference distribution 1. The vector y ∼ U(Y1 ) where Y1 = Sd . Expectation
under this distribution is denoted E1 (·).
Reference distribution 2. The vector y ∼ U(Y2 ) where
Y2 = {z ∈ Sd | z T xc = ρ̃},
p
for some −1 6 ρ̃ 6 1, and c ∈ {0, 1, . . . , N − 1}. Then y = ρ̃xc + 1 − ρ̃2 y ∗
for y ∗ uniformly distributed on a subset of Sd isomorphic to Sd−1 . Expectation
under this distribution is denoted E2 (·).
Reference distribution 1 holds true if the Yi are IID Gaussian random variables (with positive variance). In that case, the estimate p̂1 is the same as we
would get under a t-test. Reference distribution 2 is a significant narrowing of
reference distribution 1 in the direction of the ultimate reference distribution: a
point mass on y = y 0 .
3. Approximation via spherical cap volume
Here we study the approximate p-value p̂1 (ρ̂) = σd (C(y; ρ̂)). First we find the
mean squared error of this approximation over all spherical caps of the given
volume via invariance. Next we give a probabilistic interpretation which includes
the conditional unbiasedness result in Proposition 2 below. Then we give two
computational simplifications, first taking advantage of the permutation structure of our points, and then second for permutations of a binary vector. We
begin by restating the invariance principle.
Theorem 1. Let x0 , . . . , xN −1 be any points in Sd . Then
Z 1Z
N −1
N −1
2
1 X
1
1 X
kx
−
x
k
+
σ
(C(z;
t))
−
1
(x
)
dσd (z) dt
k
`
k
d
C(z;t)
N2
Cd −1 Sd
N
k,`=0
k=0
Z Z
=
kx − yk dσd (x) dσd (y)
Sd

Sd

where Cd = ωd−1 /(dωd ).
Proof. Stolarsky (1973).
Brauchart and Dick (2013) gave a simple proof of Theorem 1 using reproducing kernel Hilbert spaces. They also generalized Theorem 1 as follows.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

●

●
●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

Sd

●

●

●
●

●
●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●
●

●
●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●
●

●

●
●
●

●
●

●

●
●

●

●
●
●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●
●

●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●
●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●
●

●

●

●

●
●

●
●
●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

●
●

●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

y0

●
●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●
●

●

●

●

●

●

●

●

●

●
●

●

●

y1

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●
●
●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●

●
●

●
●

●

●

●
●

●

●

●
●

●

●

●
●

●

●
●

●

●
●

●

●

●
●

●

10

●

Fig 1: Illustration for reference distribution 1. The point y is uniformly distributed over Sd . The small open circles represent permuted vectors xk . The
point y 0 is the observed value of y. The circle around it goes through x0 and
represents a spherical cap of height y T
0 x0 . A second spherical cap of equal volume is centered at y = y 1 . We study moments of p(y, ρ̂), the fraction of xk in
the cap centered at random y.
Theorem 2. Let x0 , . . . , xN −1 be any points in Sd . Let v : [−1, 1] → (0, ∞) be
any function with an antiderivative. Then
Z 1
Z
N −1
2
1 X
1C(z;t) (xk ) dσd (z) dt
v(t)
σd (C(z; t)) −
N
−1
Sd
k=0
(3.1)
Z
Z
N
−1
X
1
= 2
Kv (xk , x` ) −
Kv (x, y) dσd (x) dσd (y)
N
Sd Sd
k,l=0

where Kv (x, y) is a reproducing kernel function defined by
Z 1
Z
Kv (x, y) =
v(t)
1C(z;t) (x)1C(z;t) (y) dσd (z) dt.
−1

(3.2)

Sd

Proof. See Theorem 5.1 in Brauchart and Dick (2013)
If we set v(t) = 1 and K(x, y) = 1 − Cd kx − yk, then we recover the original
Stolarsky formula. Note that the statement of Theorem 5.1 in Brauchart and

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

11

Dick (2013) has a sign error in their counterpart to (3.1). The corrected statement (3.1) can be verified by comparing equations (5.3) and (5.4) of Brauchart
and Dick (2013).
We would like a version of (3.1) just for one value of t such as t = ρ̂ = xT
0 y0 .
For ρ̂ ∈ [−1, 1) and  = (1 , 2 ) ∈ (0, 1)2 , let
v (t) = 2 +

1
1(ρ̂ 6 t 6 ρ̂ + 1 ).
1

(3.3)

Each v satisfies the conditions of Theorem 2 making (3.1) an identity in . We
let 2 → 0 and then 1 → 0 on both sides of (3.1) for v = v yielding Theorem 3.
Theorem 3. Let x0 , x1 , . . . , xN −1 ∈ Sd and t ∈ [−1, 1]. Then
Z

|p(y, t) − p̂1 (t)|2 dσd (y) =

Sd

N −1 N −1
1 X X
σd (C2 (xk , x` ; t)) − p̂1 (t)2 .
N2

(3.4)

k=0 `=0

Proof. See Section 11.1 of the Appendix which uses the limit argument described
above.
We now give a proposition that holds for all distributions of y ∈ Sd including
our reference distributions 1 and 2.
Proposition 1. For a random point y ∈ Sd ,
E(p(y, t)) =

N −1
1 X
Pr(y ∈ C(xk ; t)),
N

and

(3.5)

N −1
1 X
Pr(y ∈ C2 (xk , x` ; t)).
E(p(y, t) ) = 2
N

(3.6)

k=0

2

k,`=0

Proposition 2. For any x0 , . . . , xN −1 ∈ Sd and t ∈ [−1, 1], p̂1 (t) = E1 (p(y, t)).
h 1 NP
i
−1
Proof. E1 (p(y, t)) = E1
1C(y;t) (xk ) = σd (Cd (y; t)) = p̂1 (t).
N k=0
Combining Propositions 1 and 2 with Theorem 3 we find that if y ∼ U(Sd ),
as it would for IID Gaussian Yi , then p(y, ρ̂) is a random variable with mean
p̂1 (ρ̂) and variance given by (3.4) with t = ρ̂. Here ρ̂ = xT
0 y 0 is fixed while y is
random.
The right hand side of (3.4) sums O(N 2 ) terms. In a permutation analysis
1
we might have N = n! or N = m0m+m
for binary Xi , and so the computational
0
cost could be high. The symmetry in a permutation set allows us to use
Z
Sd

|p(y, t) − p̂1 (t)|2 dσd (y) =

N −1
1 X
σd (C2 (x0 , xk ; t)) − p̂1 (t)2
N
k=0

instead. This expression costs O(N ), the same as the full permutation analysis.
The cost can be reduced for binary Xi .

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

12

When the Xi are binary, then for fixed t, σd (C2 (xk , x` ; t)) just depends on
the swap distance r between xk and x` . Then
Z

|p(y, t) − p̂1 (t)|2 dσd (y) =

Sd

m
1 X
Nr V2 (u(r); t, d) − p̂1 (t)2
N r=0

for V2 (u(r); t, d) given in Lemma 1, where Nr =
pairs (xk , x` ) at swap distance r.

PN −1 PN −1
k=0

`=0

(3.7)

1(rk,` = r) counts

Theorem 4. Let x0 ∈ Sd be the centered and scaled vector from an experiment
with binary Xi of which
 m0 are negative and m1 are positive. Let x0 , dx1 , . . . , xN −1
1
be the N = m0m+m
distinct permutations of x0 . If y ∼ U(S ), then for
0
t ∈ [−1, 1], and with u(r) defined in (2.1),
E(p(y, t)) = σd (C(y 0 ; t)), and
m   
m1
1 X m0
V2 (u(r); t, d) − p̂1 (t)2 .
Var(p(y, t)) =
r
N r=0 r
Proof. There are

m0
r



m1
r



permuted points xi at swap distance r from x0 .

4. A finer approximation to the p-value
In the previous section, we studied the distribution of permutation p-values
p(y, t) with spherical cap centers y ∼ U(Sd ) and heights t = ρ̂. In this section,
we use reference distribution 2 to obtain a finer approximation to p(y 0 , ρ̂) by
studying the distribution of the p-values with centers y satisfying the constraint
hy, x0 i = hy 0 , x0 i = ρ̂. That is y has reference distribution 2, which is U(Y2 ).
Our methods also let us impose the constraint hy, xc i = hy 0 , xc i ≡ ρ̃ for
any c = 0, 1, 2, . . . , N − 1 that we like. Conditioning on hy, xc i = ρ̃ eliminates
many irrelevant y from consideration. In addition to our estimate p̂2 obtained
by xc = x0 , we consider a second choice. It is to choose xc to be the closest
permutation of x0 to y 0 . That is c = arg maxk hy 0 , xk i.
For an index c ∈ {0, 1, . . . , N − 1} conditioning as above leads to
p̃c = E2 (p(y, ρ̂)) = E1 (p(y, ρ̂) | y T xc = y T
0 xc ),

(4.1)

and our two special cases are
p̂2 ≡ p̃0 ,

and p̂3 ≡ p̃c ,

where

c = arg max hy 0 , xk i.
06k<N

(4.2)

For an illustration of reference distribution 2 see Figure 2.
Notice that p̂2 cannot go below 1/N because all of the points y that it includes
have x0 ∈ C(y; ρ̂). In fact x0 is on the boundary of this spherical cap. Since the
true value satisfies p > 1/N , having p̂2 > 1/N is a desirable property. Similarly,
p̂3 > 1/N because then xc is in general an interior point of C(y, ρ̂). We expect

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

●

●
●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●
●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●
●

●

●

●

●

Sd

●
●

●

●
●

●
●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●

●
●
●

●
●

●

●

●
●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

●

●
●

●
●

●

●

●

●

●

●
●

●

●

●

●

●
●

●
●

●

●

●

●
●

●

●
●

●

●

●

●

●

●

●

●

●
●
●

●
●

●

●
●
●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●
●

●

●

●
●

●

●
●

●

●
●

●

●

●

●

●

●

xc

●

●

●

●

●
●

●

●

●

●

●
●

●

●

●
●

●

●

●
●

●

y0

●
●

●

●

●

●
●

●

●

●

●
●

●

y2
●

●
●

●
●

●

●

●

●

●
●

●
●

●

●

●

●

●

●
●

●

y1

●

●

●

●
●

●

●
●

●

●
●

●

●
●
●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●
●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●

●
●

●

●

●

●

●

●

●

●
●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●
●

●

●

●

●

●

●
●

●
●

●

●

●
●

●

●

●
●

●

●

●
●

●

●
●

●

●
●

●

●

●
●

●

13

●

Fig 2: Illustration for reference distribution 2. The original response vector is
y 0 with y T
0 x0 = ρ̂ and x0 marked xc. We consider alternative y uniformly
distributed on the surface of C(x0 ; ρ̂) (dashed circle) with examples y 1 and
y 2 . Around each such y j there is a spherical cap of height ρ̂ that just barely
includes xc = x0 . The small open circles are permuted points xk . The fraction
of those open circles that belong to the spherical cap C(y, ρ̂) is p(y, ρ̂). We use
p̂2 = E2 (p(y, ρ̂)) and find an expression for E2 ((p̂2 − p(y, ρ̂))2 ).

that p̂3 should be more conservative than p̂2 and we see this numerically in
Section 7.
From Proposition 1, we can get our estimate p̃c and its mean squared error
by finding single and double inclusion probabilties for y.
To compute p̃c we need to sum N values Pr(y ∈ C(xk ; t) | y T xc = ρ̃) and for
p̃c to be useful we must compute it in o(N ) time. The computations are feasible
in the binary case (Xi at two levels), which we now focus on.
Let uj = xT
j = 1, 2, and let u3 = xT
1 x2 . Let the projection of y
j x0 forp
∗
2
on xc be y = ρ̃xc + 1 − ρ̃ y . Then the single and double point inclusion
probabilities under reference distribution 2 are
Z
P1 (u1 , ρ̃, ρ̂) =
1(hy, x1 i > ρ̂) dσd−1 (y ∗ ), and
(4.3)
Sd−1
Z
P2 (u1 , u2 , u3 , ρ̃, ρ̂) =
1(hy, x1 i > ρ̂)1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
(4.4)
Sd−1

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

14

where ρ̂ = hx0 , y 0 i. If two permutations of x0 are at swap distance r, then their
−1
inner product is u(r) = 1 − r(m−1
0 + m1 ) from equation (2.1).
p
Lemma 2. Let the projection of x1 onto xc be x1 = u1 xc + 1 − u21 x∗1 . Then
the single point inclusion probability from (4.3) is
(
1(ρ̃u1 > ρ̂),
u1 = ±1 or ρ̃ = ±1
P1 (u1 , ρ̃, ρ̂) =
(4.5)
σd−1 (C(x∗1 , ρ∗ )), u1 ∈ (−1, 1), ρ̃ ∈ (−1, 1)
where ρ∗ = (ρ̂ − ρ̃u1 )/

p
(1 − ρ̃2 )(1 − u21 ).

p
Proof. The projection of y onto xc is y = ρ̃xc + 1 − ρ̃2 y ∗ . Now
(
ρ̃u1 ,
u1 = ±1 or ρ̃ = ±1
p
p
hy, x1 i =
∗
∗
2
2
ρ̃u1 + 1 − ρ̃ 1 − u1 hy , x1 i , u1 ∈ (−1, 1), ρ̃ ∈ (−1, 1)
and the result easily follows.
We can now give a computable expression for p̃c and hence for p̂2 and p̂3 .
Theorem 5. For −1 6 ρ̂ 6 1 and −1 6 ρ̃ 6 1,
p̃c = E2 (p(y, ρ̂)) =

m   
1 X m0
m1
P1 (u(r), ρ̃, ρ̂)
N r=0 r
r

(4.6)

where u(r) is given in equation (2.1), P1 (u(r), ρ̃, ρ̂) is given in equation (4.5)
and ρ̃ = xT
c y0 .


Proof. There are mr0 mr1 permutations of x0 at swap distance r from xc .
From (4.6) we see that p̃c can be computed in O(m) work. The mean squared
error for p̃c is more complicated and will be more expensive. We need the double
point inclusion probabilities and then we need to count the number of pairs
xk , x` forming a given set of swap distances among xk , x` , xc .
Lemma 3. For j = 1, 2, let rj be the swap distance of xj from xc and let r3 be
the swap distance between x1 and x2 . Let u1 , u2 , u3 be the corresponding inner
products given by (2.1). If there are equalities among x1 , x2 and xc , then the
double point inclusion probability from (4.4) is

x1 = x2 = xc

1(ρ̃ > ρ̂),

1(ρ̃ > ρ̂)P (u , ρ̃, ρ̂), x = x 6= x
1 2
1
c
2
P2 (u1 , u2 , u3 , ρ̃, ρ̂) =

1(ρ̃
>
ρ̂)P
(u
,
ρ̃,
ρ̂),
x
=
x
=
6
x
1
1
2
c
1



P1 (u2 , ρ̃, ρ̂),
x1 = x2 6= xc .
If x1 , x2 and xc are three distinct points with min(u1 , u2 ) = −1, then
(
1(−ρ̃ > ρ̂)P1 (u2 , ρ̃, ρ̂), u1 = −1
P2 (u1 , u2 , u3 , ρ̃, ρ̂) =
1(−ρ̃ > ρ̂)P1 (u1 , ρ̃, ρ̂), u2 = −1.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

15

m

m

z
}|1
{ z
}|0
{
xc = ( +, +, +, +, +, · · · , +, +, +, +, −, −, −, −, · · · , −, −, −, −, − )
m

m

z
}|1
{ z
}|0
{
x1 = ( +, +, +, · · · , +, −, −, −, · · · , −, +, +, +, · · · , +, +, +, −, · · · , − )
{z
} |
{z
}
|
r1

r1

m1

m

z
}|
{ z
}|0
{
x2 = ( +, · · · , +, −, −, · · · , −, +, · · · , +, −, · · · , −, +, +, · · ·, +, −, · · · , − )
| {z }
| {z }
δ1

|

{z

r2

δ2

}

|

{z

}

r2

Fig 3: Illustration of r1 , r2 , δ1 and δ2 . The points xc , x1 and x2 each have m0
negative and m1 positive components. For j = 1, 2 the swap distance between
xj and xc is rj . There are δ1 positive components of xc where both x1 and x2
are negative, and δ2 negative components of xc where both xj are positive.

Otherwise −1 < u1 , u2 < 1, and then
P2 (u1 , u2 , u3 , ρ̃, ρ̂)


1(ρ̃u1 > ρ̂)1(ρ̃u2 > ρ̂),



R 1 ω
d−1
d−2
(1 − t2 ) 2 −1 1(t > ρ1 )1(tu∗3 > ρ2 ) dt,
=
−1 ωd−1



R 1 ωd−2

2 d−1

2 −1 1(t > ρ )σ
 −1 ω
(1
−
t
)
C x∗∗
1
d−2
2 , √
d−1

ρ̃ = ±1
ρ̃ 6= ±1, u∗3 = ±1
ρ2 −tu∗
3
1−t2

√

1−u∗2
3



dt,

ρ̃ 6= ±1, |u∗3 | < 1

where
u∗3 = p

u3 − u1 u2
p
1 − u21 1 − u22

and

ρ̂ − ρ̃uj
q
, j = 1, 2
ρj = p
1 − ρ̃2 1 − u2j

(4.7)

∗
∗
and x∗∗
2 is the residual from the projection of x2 on x1 .

Proof. See Section 11.2.
Next we consider the swap configuration among x1 , x2 and xc . Let xj be
at swap distance rj from xc , for j = 1, 2. We let δ1 be the number of positive
components of xc that are negative in both x1 and x2 . Similarly, δ2 is the
number of negative components of xc that are positive in both x1 and x2 . See
Figure 3. The swap distance between x1 and x2 is then r3 = r1 + r2 − δ1 − δ2 .
Let r = (r1 , r2 ), δ = (δ1 , δ2 ) and r = min(r1 , r2 ). We will study values of
r1 , r2 , r3 , δ1 , δ2 ranging over the following sets:
r1 , r2 ∈ R = {1, . . . , m}
δ1 ∈ D1 (r) = {max(0, r1 + r2 − m0 ), . . . , r}
δ2 ∈ D2 (r) = {max(0, r1 + r2 − m1 ), . . . , r},

and

r3 ∈ R3 (r) = {max(1, r1 + r2 − 2r), . . . , min(r1 + r2 , m, m0 + m1 − r1 − r2 )}.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

16

Whenever the lower bound for one of these sets exceeds the upper bound, we
take the set to be empty, and a sum over it to be zero. Note that while r1 = 0
is possible, it corresponds to x1 = xc and we will handle that case specially,
excluding it from R.
The number of pairs (x` , xk ) with a fixed r and δ is
  




m0
m1
m0 − δ1
m1 − δ2
m0 − r1
m1 − r1
c(r, δ) =
. (4.8)
δ1
δ2
r1 − δ1
r1 − δ2
r2 − δ1
r2 − δ2
Then the number of configurations given r1 , r2 and r3 is
X X
c(r, δ)1(r3 = r1 + r2 − δ1 − δ2 ).
c(r1 , r2 , r3 ) =

(4.9)

δ1 ∈D1 δ2 ∈D2

We can now get an expression for the expected mean squared error under
reference distribution 2 which combined with Theorem 5 for the mean provides
an expression for the mean squared error of p̃c .
Theorem 6. For −1 6 ρ̂ 6 1 and −1 6 ρ̃ 6 1,

m   
X
m0
m1
1
P2 (1, u(r), u(r), ρ̃, ρ̂)
1(ρ̃
>
ρ̂)
+
2
r
r
N2
r=1
m   
X
m0
m1
(4.10)
+
P1 (u(r), ρ̃, ρ̂)
r
r
r=1

X X X
+
c(r1 , r2 , r3 )P2 (u1 , u2 , u3 , ρ̃, ρ̂)

E2 (p(y, ρ̂)2 ) =

r1 ∈R r2 ∈R r3 ∈R3 (r)

where P2 (·) is the double inclusion probability in (4.4) and c(r1 , r2 , r3 ) is the
configuration count in (4.9).
Proof. See Section 11.3 of the Appendix.
In our experience, the cost of computing E2 (p(y, ρ̂)2 ) under reference distribution 2 is dominated by the cost of the O(m3 ) integrals required to get the P2 (·)
values in (4.10). The cost also includes an O(m4 ) component because c(r1 , r2 , r3 )
is also a sum of O(m) terms, but it did not dominate the computation at the
sample sizes we looked at (up to several hundred).
5. Generalized Stolarsky Invariance
Here we obtain the results for reference distribution 2 in a different way, by
extending the work by Brauchart and Dick (2013). They introduced a weight
on the height t of the spherical cap in the average. We now apply a weight
function to the inner product hz, xc i between the center z of the spherical cap
and a special point xc .

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

17

Theorem 7. Let x0 , . . . , xN −1 be arbitrary points in Sd and v(·) and h(·) be
positive functions in L2 ([−1, 1]). Then for any x0 ∈ Sd , the following equation
holds,
N −1
2
1 X
v(t)
1C(z;t) (xk ) dσd (z) dt
h(hz, x i) σd (C(z; t)) −
N
−1
Sd
k=0
Z Z
N
−1
X
1
= 2
Kv,h,x0 (xk , x` ) +
Kv,h,x0 (x, y) dσd (x) dσd (y)
N
Sd Sd

Z

1

Z

0

(5.1)

k,`=0

−

N −1 Z
2 X
Kv,h,x0 (x, xk ) dσd (x)
N
Sd
k=0

where Kv,h,x0 : Sd × Sd → R is a reproducing kernel defined by
Z 1
Z
Kv,h,x0 (x, y) =
v(t)
h(hz, x0 i)1C(z;t) (x)1C(z;t) (y) dσd (z) dt.

(5.2)

Sd

−1

Proof. See Section 11.4 of the Appendix.
Remark. We will use this result for x0 = xc , where xc is one of the N given
points. The theorem holds for general x0 ∈ Sd , but the result is computationally
and statistically more attractive when x0 = xc .
We now show that the second moment in Theorem 6 holds as a special limiting
case of Theorem 7. In addition to v from Section 3 we introduce η = (η1 , η2 ) ∈
(0, 1)2 and
hη (s) = η2 +

η1 ( ωωd−1
(1
d

1
1(ρ̃ 6 s 6 ρ̃ + η1 )
− s2 )d/2−1 )

(5.3)

Using these results we can now establish the following theorem, which provides the second moment of p(y, ρ̂) under reference distribution 2.
Theorem 8. Let x0 ∈ Sd be the centered and scaled vector from an experiment
with binary Xi of which
 m0 are negative and m1 are positive. Let x0 , x1 , . . . , xN −1
1
be the N = m0m+m
distinct permutations of x0 . Let xc be one of the xk and
0
define p̃c by (4.1). Then
E2 (p̃c (y, ρ̂)2 ) =

N −1 Z
1 X
1(hy, xk i > ρ̂)1(hy, x` i > ρ̂) dσd−1 (y ∗ )
N2
Sd−1
k,`=0

where y = ρ̃xc +

p

1 − ρ̃2 y ∗ .

Proof. The proof uses Theorem 7 with a sequence of h defined in (5.3) and v
defined in (3.3). See Section 11.5 of the appendix.
This result shows that we can use the invariance principle to derive the second
moment of p(y, ρ̂) under reference distribution 2. The mean square in Theorem 8
is consistent with the second moment equation (3.6) in Proposition 1.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

18

6. Two-sided p-values
In statistical applications it is more usual to report two-sided p-values. A conservative approach is to use 2 min(p, 1 − p) where p is a one-sided p-value. A
sharper choice is
N −1
1 X
p=
1(|xT
k y 0 | > |ρ̂|).
N
k=0

This choice changes our estimate under reference distribution 2. It also changes
the second moment of our estimate p̂1 .
The two-sided version of the estimate p̂1 (ρ̂) is 2σd (C(y; |ρ̂|)), the same as if
we had doubled a one-sided estimate. Also E1 (p) = p̂1 in the two-sided case.
We now consider the mean square for the two-sided estimate under reference
distribution 1. For x1 , x2 ∈ Sd with u = xT
1 x2 , the two-sided double inclusion
probability under reference distribution 1 is
Z
Ṽ2 (u; t, d) =
1(|z T x1 | > |t|)1(|z T x2 | > |t|) dσd (z).
Sd
T

Writing 1(|z xi | > |t|) = 1(z T xk > |t|) + 1(z T (−xk ) > |t|) for k = 1, 2 and
expanding the product, we get
Ṽ2 (u; t, d) = 2V2 (u; |t|, d) + 2V2 (−u; |t|, d).
By replacing V2 (u, t, d) with Ṽ2 (u, t, d) and p̂1 (t) with 2σd (C(y; |t|)) in Theorem 4, we get the variance of two-sided p-values under reference distribution 1.
To obtain corresponding formulas under reference distribution 2, we use the
T
usual notation. Let uj = xT
j x0 for j = 1, 2, and let u3 = x1 x2 . Let the projection
p
of y on xc be y = ρ̃xc + 1 − ρ̃2 y ∗ . Now
Z

P̃1 (u1 , ρ̃, ρ̂) =
1 |hy, x1 i| > |ρ̂| dσd−1 (y ∗ ), and,
(6.1)
d−1
ZS


P̃2 (u1 , u2 , u3 , ρ̃, ρ̂) =
1 |hy, x1 i| > |ρ̂| 1 |hy, x2 i| > |ρ̂| dσd−1 (y ∗ ) (6.2)
Sd−1

are the appropriate single and double inclusion probabilities.
After writing 1(|hy, xk i| > |ρ̂|) = 1(hy, xk i > |ρ̂|) + 1(hy, −xk i > |ρ̂|) for
k = 1, 2 and expanding the product, we get
P̃1 (u1 , ρ̃, ρ̂) = P1 (u1 , ρ̃, |ρ̂|) + P1 (−u1 , ρ̃, |ρ̂|),

and

P̃2 (u1 , u2 , u3 , ρ̃, ρ̂) = P2 (u1 , u2 , u3 , ρ̃, |ρ̂|) + P2 (−u1 , u2 , −u3 , ρ̃, |ρ̂|)
+ P2 (u1 , −u2 , −u3 , ρ̃, |ρ̂|) + P2 (−u1 , −u2 , u3 , ρ̃, |ρ̂|).
Changing P1 (u1 , ρ̃, ρ̂) and P2 (u1 , u2 , u3 , ρ̃, ρ̂) to P̃1 (u1 , ρ̃, ρ̂) and P̃2 (u1 , u2 , u3 , ρ̃, ρ̂)
respectively in Theorems 5 and 6, we get the first and second moments for twosided p-values under reference distribution 2.
For a two-sided p-value, p̂3 is calculated with xc where c̃ = arg maxk |hy 0 , xk i|.
For m0 = m1 , c̃ = c = arg maxk hy 0 , xk i, but the result may differ significantly
for unequal sample sizes.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

19

7. Numerical Results
We consider two-sided p-values in this section. The main finding is that the root
mean squared error (RMSE) of p̂2 under reference distribution 2 is usually just
a small multiple of p̂2 itself.
First we evaluate the accuracy of p̂1 , the simple spherical cap volume approximate p value. We considered m0 = m1 in a range of values from 5 to 200. The
values p̂1 ranged from just below 1 to 2 × 10−30 . We judge the accuracy of this
estimate by its RMSE. Under distribution 1 this is (E(p̂1 (ρ) − p(y, ρ))2 )1/2 for
y ∼ U(Sd ). Figure 4a shows this RMSE decreasing towards 0 as p̂1 goes to 0
with ρ going to 1. The RMSE also decreases with increasing sample size, as we
would expect from the central limit theorem.
As seen in Figures 4a and 4b, the RMSE is not monotone in p̂1 . Right at
p̂1 = 1 we know that RMSE = 0 and around 0.1 there is a dip. The practically
interesting values of p̂1 are much smaller than 0.1, and the RMSE is monotone
for them.
A problem with p̂1 is that it can approach 0 even though p > 1/N must hold.
The distribution 1 RMSE does not reflect this problem. By studying E2 ((p̂1 (ρ)−
p(y, ρ))2 )1/2 , we get a different result. In Figure 4c, the RMSE of p̂1 under
distribution 2 reaches a plateau as p̂1 goes to 0.
The estimator p̂2 = p̃0 performs better than p̂1 because it makes more use
of the data, and it is never below 1/N . As seen in Figure 4d, the RMSE of
p̂2 very closely matches p̂2 itself as p̂2 decreases to zero. That is, the relative
error |p̂2 − p|/p̂2 is well behaved for small p-values. In rare event estimation,
that property is known as strong efficiency (Blanchet and Glynn, 2008) and can
be very hard to achieve. Here as p̂2 decreases to the granularity limit 1/N , its
RMSE actually decreases to 0. Eventually the distance from y 0 to x0 is below
the minimum interpoint distance among the xk and then, for a one-sided test,
p̂2 = p = 1/N .
The estimators p̂1 and p̂2 , do not differ much for larger p-values as seen in
Figure 5a. But in the limit as ρ̂ → 1 we see that p̂1 → 0, while p̂2 approaches
the granularity limit 1/N instead.
Figure 5b compares the RMSE of the two estimators under distribution 2. As
expected, p̂2 is more accurate. It also shows that the biggest differences occur
only when p̂1 goes below 1/N .
To examine the behavior of p̂2 more closely, we plot its coefficient of variation in Figure 6. We see that the relative uncertainty in p̂2 is not extremely
large. Even when the estimated p-values are as small as 10−30 the coefficient of
variation is below 5.
In Section 4, we mentioned another choice for xc . It was p̂3 = p̃c , where xc
is the closest permutation of x0 to y 0 . Figure 2.7 in He (2016) compares p̂3 to
p̂2 in some simulations. As expected, p̂3 tends to be larger (more conservative)
than p̂2 , though it does sometimes come out smaller. Figure 2.8 of He (2016)
compares the RMSE of p̂3 to p̂2 . The upward bias of p̂3 gave it a much larger
RMSE.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

−20

−30

−30

−20

−10

●●●●●
●●●●●●●●●● ●●●●●●●
●●●●●●●
●●
●●●●●●●
●
●●●●●●
●
●●●●●
●●●●●
●
●●●●●
●●●●●
●
●●●●●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●●●●●●●●●●●●●●●
●●●●●●●●●●
●
●●●●●●●●●●
●●●●●●●●●●●
●●●●●●●●●●●
●●●●●●●●●●
●●●●●●●●●
●●●●●●●●
●●●●●●●
●●●●●●●●●●●●●●●●●●●●
●●●
●●●●●
●●●●
●●
●
●
●
●●
●●●
●●
●●●
●●●
●
●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●●●
●
●●
●●●●●
●●●
●●
●
●●●●●
●●
●●
●●●●
●●
●
●●●
●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●●●
●●
●●
●●●●●●●●●●●●●●●●●
●●●●●●●
●●●
●●●●●●●●●●●●●●●●●●●●
●●
●●●●●●●●●●●
●●●●●
● ●
●●●
●●●●●●
●●●●●●●●●
●●●
●●●●
●●
●
●●
●●●●
●●●●●●●●
●●●●
●●
●●
●●
●●●●
●● ● ●
●●●●
●●
●●
●●●
●●
●● ●
●●●●
●●●
●●
●●●●●●●●●●●●
● ●
●●●
●●
●
●
●●●
●●
●●●●●●●●●●●●●
● ●
●●●●●●●●●●●●●●●●●●●●●●●●●●
●●
●●
●●●●●●●●●
●●●●
●●●●●●●●●●●●●●●●●●
●●●●●●●
●●
●
●●●●●
●●●●●●●●●●●
●●●●●
●● ●●●
●●●
●●●●
●●
●●●●●●●●●
●●●●
●●
●● ●●
●●●●
●●●
●●●
●● ● ●
●●●
●● ●●●
●
●●●
●●
● ●●
●●●●●●●●●●●●●●●●●●●●●●●●
●●●
●● ●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●●
●●●●●
●
●●●
● ●●●
●●●●●●●●●●●●●
●●●●●●
●●●
●●●●
●●
● ●
●
●
●
●
●
●●●●●●●●●●
●
●●●●
●
●●
●● ●
●●
●
● ●●
●●●●●●●●●●●●●●●
●●●
●●
●●●
● ● ●
●●
●●●●●●●●
●●●
●●●
●● ●
●●●
●●
●●
●●●●●
●●● ●● ● ●
●●●
●●●
●●●●
●●● ●●● ●●
●●
●● ●
●●
●
●●●●
●●●●
●●●
●●
●● ● ●
●●
●●●
●
●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●●
●●
●● ●
●●●
●●
●
●●●●●●●●●●●●●●●●●
●●●●●●●
●●
● ●●● ●●● ●●●●●
●●●●●●●●●●●
●●●●●
● ● ●●
●
●
●
●●
●●●●
●●
● ●● ●●● ●●●●
●●●
●● ●
●●
●
●
●
●
●
●●
● ●
● ● ●● ●●●
●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●●
●●●●●●●●●●●●
●●●●●
●●
●●
●● ●●● ●● ●● ●●●
●●●●●●●●●●●●●●●●
●●●●
●●
●●
●●●●●●●●
●●●
● ● ●
●●●
● ●●
●●●●●
●●●
●●
●●●
● ●● ●●
●●●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●
●●
●●
●●●●●●●●●●●●
●●●●
●
●●●●●
● ● ●●
●
●●
●●
●●●
●
●●●●
● ●● ●●
●●
●●
●●●
●●●
●
●●
●●●
●●●
● ● ●●
●
●
●●
●●
●●●
●●●
●
●● ●●● ●●
●●
●
●●
●●
● ●●●
●●
●● ●● ●
● ●
●
●●
●● ● ● ●●●●●
● ●
●●
●●●●●●●●●●●●●●●●●●●●●●
● ● ●●●● ●●
●●
●●●
●●●●●●
●
● ●
●●
●●●●●
●●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●
●●
●
●●●●
●●
●●●●●●●●●●●●●
●●●●●●
●●
●●●
●●
● ●
●
●●●●
●●●
●●
●●
● ●
●●
●●●
●●●
●
●
●●●
●
●
● ● ●
●
●
●●●
●●●
●
● ●●
● ● ●
●●
●●
●
●
●
●●
●●
● ●●● ●●
●●
●
●●
●●
●
●
●●
●
●
●
●
●
● ●
●●
●
●
●
●●
●
●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●
●
●●●●●●
●●●●●●●●●●●●●●
●●
●
●
●●● ●
●●●●
●●●●
●●
●
●
●
●●●
●
●●
●
●●●
●
●●
●
●
●
●●
●
●
●●
●
● ●
●
●●
●
●
● ● ●
●
●●
●
●
●
●
●●
●
●
● ● ●
● ●
●●●
● ●
●
●

−2

(m0, m1)
●

●

●

●

●

●

●

●

●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
100,100
200,200

log10(RMSE1)

log10(RMSE1)

−10

●●●
●●
●●●●●
●●●
●●
●●●
●●
●●
●●●
● ●●
●●●
●●● ●●●●●●●●●●●
●●
●●●●
●● ●●●●●● ●
●●
●
●●
●●
●●●●●
●●●●● ●●
●●
●●
●●
●●●●●
●●●●●●●●●●●
●●● ●
●●
●●●●●●
●●●●●
●●
●● ●●●●●●
●●●●
●●
●●●●●
●● ●●●●
●● ●
●●
●●●●●●●●
●●●●
●
●●●●●
●●
●●●●●●
●●●●
●●●
●●
●●●●●
●
●
●
●
●
●
●
●
●
●●● ●●● ●●●● ●
●●●●●●●
●●
●●●●
●
●●
●
●
●●
●●●●●●
●●●●● ●● ●●●● ●
●●
●●●●●●●●
●●● ●● ●●● ●●● ●
●●
●●●●●●●●●● ●●●
●●●●●
●●
●
●
●
●● ● ●
●
●●●
●● ●●●●●
●●
●
● ●● ●●
●● ●●●●●●
●●
●● ●●●●●●
●●● ●● ●●
●●
●● ●●
●●
●●●●●● ●● ●●
●
●● ●●●
●●
●
●
●
●
●
●
●
●
●
●
●● ●● ●
●●●
●●
●●
●
●●●●●●
●●
●●
●●●●●●●●●
●
●●●●●●
●●
●●
●●
●●
●●●●●●●●
●
●●●●●●●
●●
●●
●
●
●●●●●●●●
●●
●●
●●●●●●●●
●●
●●
●●●●●●
●
●●
●●
●
●●●●●●●
●●●●●●
●●
●●
●
●●
●●
●●
●●●●●●●●
●●
●●
●● ●●●●
●● ●●●●
●●
●●
●●
●●
●●●●●●●
●● ●●●●
●●
●●
●
●●
●● ●●●●●
●●
●●
●● ●●●
●●●●●●●
●●
●●
●● ●●●●
●●
●●
●●●●●●
●●
●●
●● ●●●●●
●●
●●
●● ●●●
● ●●
●●
●●
●●●● ●
●● ●●●●●
●●
●●
●●●●●● ●●
●●
●●
●●
●
●● ●●●●●●●
●●
●●
●●
●● ●●●●
●●
●●
●
●
●
●
●
●
●
● ●●●●●●● ●
●
●
●●
●
●●
●●
●● ●●
●
●●
●●
●●
●●●●●●
●●
●● ●●●●
●●
●●
●●
●●
●●
●● ●●●●
●● ●●●●
●●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●●● ● ●
●●
●●●
●●
●●
●● ●●●●
●●
●●
●●
●● ●●
● ●●
●●
●●
●●
●●
●●●●●●
●● ●●●
●●
●●
●●
●●
●●
●●
●● ●●●●●
●●
●●
●●
●● ●●●●
●●
●●
●●
●●
●●
●●
●●
●● ●●●●●
●●
●●
●●
●● ●●●●
●●
●●
●●
●●
●●
●●
●●
●●
●●●●
●●
●●
●●
●●
●●
●●
●●●●● ●●●●
●●
●●
●●
●●
●●
●●
●●
●●●●●●●●●●
●●
●●
●●
●●
●●
●●
●●
●●●●●●●●●●
●●
●●
●●
●●
●●●
●●
●●
●●●●
●●
●●
●●
●●
●●
●●
●●
●●●●
●●
●●
●●
●●
●● ● ●
●●
●●
●●
●●
●●●●● ●●●●
●●
●●
●●
●●
●●
●●●●●● ● ●
●●
●●
●●
●●
●●
●●
●●
●●●●●● ●●●●
●●
●●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●● ●●● ●●●●
●●
●●
●●
●●
●●
●●
●●
●● ●●●● ● ●
●●
●●
●●
●●
●●
●● ●●●● ●●●●
●●
●●
●●
●●
●●
●●
●●
●● ●● ●●●●
●●
●●
●●
●●
●●
●● ●●●● ● ●
●●
●●
●●
●●
●●
●●
●●
●● ●●●● ●●●●
●●
●●
●●
●●
●● ●●● ● ●
●●
●●
●●
●●
●●
●●
●● ●●● ●●●●
●●
●●
●●
●● ●●●● ● ●
●●
●●
●●
●●
●●
●●
●● ●●● ●●●●
●●
●●
●●
●● ●●● ● ●
●●
●●
●●
●●
●●
●● ●●● ●●●●
●●
●●
●●
●●
●● ●●●● ●●●●
●●
●●
●● ●●
●●
●●
●●
●●
●● ●●●●● ●●●●
●●
●●
●
●●
●●
●●
●●
●●●● ●●●●
●●
●●
●●
●
●
●
●
●
● ●● ● ●●
●●
●●
●●
●●
●●
●●●● ●● ●
●●
●●
●●
●●
●●
●● ● ● ●●
●●
●●
●● ●● ●●●●
●●
●●
●●
●● ●● ●● ●
●●
●●
●●
●● ●● ● ●●
●●
●●
●●
●● ●● ●●●●
●●
●●
●● ●● ●● ●
●●
●●
●●
●●
●● ● ● ●●
●●
●● ●● ●● ●
●●
●●
●● ●● ● ●●
●●
●● ●● ●●●●
●●
●●
●● ● ●● ●
●●
●● ●● ● ●●
●●
●● ●● ●● ●
●●
●●
●● ●● ●● ●●
●●
●
●
●
●
●
● ●
●●
●● ●● ●● ●●
●●
●● ●●
● ●●
●●
●●
●● ●●
●● ●
●●
●
●●
● ●●
●●
●● ●●●
●●
●●●●
●●
●●
●● ●
●●
● ●●
●●
●●
●●●●
●●
●●
●●
●●
●●●●
●●
●●
●●●●
●●
●●
●●
●● ●
●●
● ●●
●●
●●
●●
●●●●
●●
●●
●●●●
●●
●●
●
●●
●
● ●●
●●
●●●●
●●
●● ●
●●
●● ●●
●●
●●
●● ●●
●●
●● ●●
●● ●●
●● ●●
●● ●●
●● ●●
●● ●●
●● ●●
●● ●●
●
●
●
●

−4

−6

●●●●●●●●

●

●

●

●

●

●

●

●

●
●

−8
−2.0

−1.5

−1.0

−0.5

0.0

log10(phat1)

(b) RMSE1 (p̂1 ) zoomed.
0

(m0, m1)
●
●
●
●

●

●
●
●

●

●

●

●

●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
70,70
100,100

log10(RMSE2)

log10(RMSE1UNDER2)

−40

●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
100,100
200,200

●

●
●

●

●

●
●
●

●

●

0

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●
●●
●●●
●●
●● ●●
●
●●●●●●●●●● ●●●
●
●●
● ●
● ●●
●
● ●●
●● ●
●●●●●●●●●
● ●●
●
●
● ●
●
●
●●
●
●●●●●●●●●
●
●
●

●

●

●

(a) RMSE1 (p̂1 ).

−20

(m0, m1)

●

log10(phat1)

0

20

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
● ●
●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●●
●●
●
● ●
●● ●
●● ●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●● ●
●
●● ●
●
●●
●
●
●●
●
●
●
●
●
●●
●
●
●
● ●
●
● ●
●●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●● ●

−20

●
●
●
●
●

●
●●

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
● ●
●

−40

(m0, m1)

●
●
●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
70,70
100,100

●
●

●

●

●

●

●

●

●
●

●
●
●
● ● ● ● ● ● ●

●

●

●

−60

−60

●
●
●
●
●
●
●
●
●
●
●
●

−120

−90

−60

−30

0

−60

−40

log10(phat1)

(c) RMSE2 (p̂1 ).

−20

0

log10(phat2)

(d) RMSE2 (p̂2 ).

Fig 4: RMSEs for p̂1 and p̂2 under reference distributions 1 and 2. The x-axis
shows the estimate p̂ as ρ varies from 1 to 0. Here m0 = m1 . Plots with m0 6= m1
are similar.

8. Comparison to saddlepoint approximation
The small relative error property of p̂2 is similar to the relative error property
in saddlepoint approximations. Reid (1988) surveys saddlepoint approximations
and Robinson (1982) develops them for permutation tests of the linear statistics
we have considered here. When the true p-value is p, the saddlepoint approximation p̂s satisfies p̂s = p(1 + O(1/n)). Because we do not know the implied
constant in O(1/n) or the n at which it takes effect, the saddlepoint approximation does not provide a computable upper bound for the true permutation
p-value p.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations
1.00

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●●
●●
●●
●● ●●
●●●●●●●●●●● ●●●
●
●●
●
●●
● ●
● ●●
●●
●●●●●●●●●● ●
● ●●
●
●
● ●
●
●
●
●
●●●●●●●●●●

−20

●

●

(m0, m1)
●
●
●
●
●

●
●

●

−40

●●●●●●●●●

●

●

●

●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
70,70
100,100

●

●
●

●
●

●

RMSE2/RMSE1UNDER2

log10(phat2)

0

●

●

●

●

●
●

0.75

●
●
●
●
●
●
●
●
●
●
●●
●●
●
●●●
●
●
●
●●
● ● ● ● ● ● ● ●●
●
●●
●● ● ●
●
● ●●
●
●
●●●
●
●
● ●● ●●
●●●●
●●
●●
●●
●
●
●
●●●●
●●
●
● ● ●
●●
●
●
●
●
● ● ●●
●
●
●
●●● ● ●
●
● ● ●●
● ●●
●● ● ● ●
● ●
●
●● ●●● ●
●●
●● ● ●
●
●● ●●
●
● ●
●
●
●
●●
●
●
●
●●
●
● ●
●
● ●
●

● ●

0.50

●

●
●
●

●

●
●

0.25

●
●

●

●

●

●

0.00

●

●

−120

●
●

●

−90

log10(phat1)

●

●

●

●
●

●
●

●

●

●

5,5
10,10
15,15
20,20
30,30
40,40
50,50
70,70
100,100

●

●

0

●

●

●

●

−30

●

●

●

−60

●
●

● ●

●

●

−90

●
● ●

●

−120

●
●

●

●

●

●

●
●●

●

●
● ● ● ● ● ● ● ●

●●

●
●

●

−60

(m0, m1)
●

● ●

●
●

●

21

●
●●

●

●
●
●●

−60

●
●
●●

−30

●

●

●

●

● ● ●
● ●
● ●
● ●
●
● ●
● ●
●
●
●
●
●
●
●
●
● ●

0

log10(phat1)

(a)

(b)

Fig 5: Comparison of p̂1 and p̂2 . In (a), log10 (p̂2 ) is plotted against log10 (p̂1 ) for
varying ρ’s. The black line is the 45 degree line. In (b), the ratio of RMSEs for
p̂1 and p̂2 is plotted against log10 (p̂1 ). The x-axis is log10 (p̂1 ).
●

●

●
●

●
●

RMSE2/phat2

10

(m0, m1)

●
●

●
●
●
●

●
●

●
●●
●
●
●

5

●
●
●
●
●
●
●

● ●
●

●
●●
● ●
●
●
●

●
●

●
●
●

●

0

●
●
●
●

−60

●

●

●
●

●
●
●
●
●

●
●
●
●
●
●

−40

5,5
10,10
15,15
20,20
30,30
40,40
50,50
70,70
100,100

●

●●●
●
●
●
●
● ● ● ●
●●●●
● ●
●
●
●● ● ●
●
●●
● ●● ● ●●●
●
● ●● ●
●●
● ●
● ●●●
●●● ●
●
●
●
●
●●●●
●●
●
●
●
●
●● ● ●●
●●
●
●
●
● ●
●
●●●
●●
●
●
●
●
●
●
● ●
●
●
●●
●
●
●
●
●●●
●
●
●
●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
● ●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●●
●●●●
●
●
●
●
●
●
●
●
●

−20

0

log10(phat2)

Fig 6: The coefficient of variation for p̂2 with varying ρ’s.

Figure 7 compares our estimates to each other and those of the saddlepoint
approximation, equation (1) from Robinson (1982). The simulated data have the
Exp(1) distribution under the control condition and the 2 + Exp(1) distribution
under the affected condition. The sample sizes were m0 = m1 = 10 making it
feasible to compute the exact permutation p-value for hundreds of examples.
In each case we ran 500 independent simulations. Cases with perfect separation
were excluded: the saddlepoint approximation is numerically unstable then, and

He, Basu, Zhao & Owen/Stolarsky invariance for permutations
phat1

phat2

phat3

saddle.point

22

60

40

20

count

0

60

40

20

0
0.1

10.0

0.1

10.0

phat/p.exact

iid

Fig 7: Simulation results p̂/p as described in the text, for Y0,i ∼ Exp(1), and
iid

Y1,i ∼ Exp(1) + 2.

one can easily detect that the minimum Y value in one group is larger than
the maximum in the other group, showing that p = 1/N . In every instance we
compared two-sided p-values. Chapter 2 of He (2016) considers simulations from
some other distributions. The control condition data are t(5) , N (0, 1) and U(0, 1)
while the affected condition data are shifted versions of these distributions.
In these simulations, the naive spherical cap estimator p̂1 , with no good
relative error properties, is consistently least accurate and is often much smaller
than the true p. The saddlepoint estimate is very accurate but tends to come
out smaller than the true p. The estimators p̂2 and p̂3 are less likely to be
below p than the saddlepoint estimate, and by construction, they are never
below the granularity limit. Qualitatively similar results happened for all of the
distributions. The accuracy of all of these p-value estimates tends to be better
for ligher tailed Yi .
We can also construct Z scores, Z2 = (p − p̂2 )/RMSE2 and a similar Z3 . If
these take large values, then it means that p̂ is too small and, moreover, that
our computed RMSE does not diagnose it. The largest Z scores we observed are
.
in Table 1. The largest Z values arose for exponential data with p = 0.89 and

He, Basu, Zhao & Owen/Stolarsky invariance for permutations
Dist’n Y0,i
t(5)
Exp(1)
U(0, 1)
N (0, 1)

max Z2

max{Z2 | p̂2 < 0.1}

max Z3

23

max{Z3 | p̂3 < 0.1}

26.7
7.55
3.49
3.07

1.91
31.5
3.87
7.55
7.76
7.76
2.45
5.87
2.61
2.78
3.07
2.78
Table 1
Maximal Z scores observed for p̂2 and p̂3 in 500 independent replications.
First author

m1

m0

N =

m1 +m0 
m1

11
18
3.5 × 107
29
14
7.9 × 1010
50
22
1.8 × 1018
Table 2
Sample sizes for three microarray studies.

Zhang
Moran
Scherzer

.
.
p̂2 = 0.78 = p̂3 . Such large p-values are not very important and so maximal Z
scores are also shown among estimated p-values below 0.1.
What we find is that the Z values are not very extreme. This suggests that it
might be feasible to get a conservative p-value estimate by adding some multiple
of the distribution 2 RMSE to p̂2 .
9. Data comparisons
Three data sets on Parkinson’s disease were used by Larson and Owen (2015)
and investigated in Chapter 6 of He (2016). They come from Scherzer et al.
(2007), Moran et al. (2006) and Zhang et al. (2005). Table 2 shows their sample
sizes.
For this comparison, there were 6180 gene sets from v5.1 of mSigDB’s gene
set collections. Curated gene sets and Gene Ontology gene sets were used. The
gene sets ranged in size from 5 to 2131 genes with an average size of 93.08 genes.
Slightly different versions of the gene sets were used in Larson and Owen (2015).
Ground truth estimates of two-sided p values for linear test statistics were
obtained using M = 106 permutations. When the estimate was below 10−4 , the
number M was increased to 107 . The Zhang data set had the smallest sample
size and had no gene sets significant at below 0.01 and so we do not compare
estimates for this gene set.
Table 3 gives correlations between log10 estimated and gold-standard p-values
for these genes. From Table 3, we see that p̂1 , p̂2 and p̂3 have nearly the same
correlations with the gold standard; indeed they correlate highly with each other.
They correlate with the gold standard estimate much more closely than the
saddlepoint estimator does. Figures in Chapter 6 of He (2016) give scatterplots.
These show the saddlepoint estimator is biased low and p̂3 is biased slightly
high. Statistics p̂1 and p̂2 are quite close, possibly because none of the gene sets
has a very small p-value.
While saddlepoint methods have a very desirable relative error property they
do have some numerical issues. Table 4 shows some timing data. We also got in-

He, Basu, Zhao & Owen/Stolarsky invariance for permutations
Data

Correlation

p̂1

p̂2

p̂3

p̂saddle

Moran
Moran

Pearson
Kendall

0.9997
0.9856

0.9997
0.9856

0.9997
0.9866

0.9937
0.9395

Scherzer
Scherzer

Pearson
Kendall

0.9997
0.9869

0.9997
0.9869

0.9997
0.9870

0.9837
0.8937

Moran Low
Moran Low

Pearson
Kendall

0.9501
0.8283

0.9504
0.8283

0.9653
0.8578

0.7125
0.5411

Scherzer Low
Scherzer Low

Pearson
Kendall

24

0.9940
0.9940
0.9947
0.8652
0.9429
0.9429
0.9429
0.7714
Table 3
The table gives both Pearson and Kendall correlations over 6180 gene sets, between
estimated log10 (p)-values and gold standard log10 (p)-values for both the Moran and
Scherzer data sets. The designation ‘Low’ refers to only the 190 gene sets with p-values in
the interval (10−5 , 10−4 ) for Moran, or the 15 gene sets with p-values in the interval
(10−5 , 10−3 ) for Scherzer.
Data Set

Saddle

Zhang
Moran
Scherzer

0.0631
0.0894
0.1394

p̂1

p̂2

0.0024
0.0031
0.0029
0.0037
0.0034
0.0045
Table 4
Average, over 6180 gene sets of the running time in seconds
p̂j , j = 1, 2, 3. Total time ranges from under 1/4 minute

p̂3
0.0032
0.0038
0.0047
for saddlepoint p-values and
to just over 14 minutes.

finite values for 95 of the gene sets on the Zhang data. It might be a convergence
issue or possibly a flaw in how we implemented saddlepoints.
10. Discussion
We have constructed approximations to the permutation p-value using probability and geometry derived from discrepancy theory. A rigorous upper bound
for p could be attained using L∞ spherical cap discrepancies instead of the L2
version, but computing such discrepancies is a major challenge. Narcowich et al.
(2010) give upper bounds for the L∞ spherical cap discrepancy, in terms of averages of a great many harmonic functions at the points xi . For our application
we need bounds for spherical caps of a fixed volume (under distribution 1) and
of fixed volume and constrained location (under distribution 2) and those go
beyond what is in Narcowich et al. (2010).
Many other approximation methods have been proposed for permutation
tests. For instance, Zhou et al. (2009) fit approximations by moments in the
Pearson family. Larson and Owen (2015) fit Gaussian and beta approximations
to linear statistics and gamma approximations to quadratic statistics for gene
set testing problems. Knijnenburg et al. (2009) fit generalized extreme value
distributions to the tails of sampled permutation values.
None of these approximations come with an all inclusive p-value that accounts for both numerical uncertainty of the estimation and sampling uncer-

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

25

tainty behind the original data. IID sampling of permutations does come with
such a p-value if we add one to numerator and denominator as Barnard (1963)
suggests. Then the Monte Carlo p-value estimate p̂ is actually a conservative
p-value in it’s own right: Pr(p̂ 6 u) 6 u for 0 6 u 6 1 under the null hypothesis.
However, that method cannot attain very small p-values, and so a gap remains.
We have employed reference distributions in an effort to address this gap. We
select a set Y containing y 0 and find the first two moments of p(y, ρ̂) for y ∼
U(Y). If the data y 0 were actually sampled from our reference distribution, then
we could get an all inclusive conservative p-value via the Chebychev inequality.
To illustrate the Chebychev inequality, let µ = E(p(y, ρ̂)) and σ 2 = Var(p(y, ρ̂))
for the observed value ρ̂ = xT
0 y 0 and for random y ∼ U(Y) for some reference
set Y. Then Pr(p > µ + λσ) 6 1/(1 + λ2 ) for any λ > 0. Under this model,
p∗ = µ + λσ + 1/(1 + λ2 ) is a conservative p-value. Minimizing p∗ over λ reduces to solving 2λ = σ(1 + λ2 )2 . For small p we anticipate λ  1 and hence
λ0 = (2/σ)1/3 will be almost as good as the optimal λ we could find numerically.
That choice leads to p∗ 6 µ + (21/3 + 2−2/3 )σ 2/3 .
For a numerical illustration, consider µ = 10−30 and σ = 3 × 10−30 , roughly
describing the small p-value estimates from the case m0 = m1 = 70. Then
p∗ 6 4 × 10−20 , much larger than µ and yet still very small, likely small enough
to be significant after multiplicity adjustments.
The reference distributions describe a set in which p(y 0 , ρ̂) is known to lie.
Reference distribution 1 applies for Gaussian Yi and reference distribution 2 is
for Gaussian Yi after conditioning on y T x0 = y T
0 x0 . Of course, the data will
not ordinarily be exactly Gaussian. The numerical illustration uses a Chebychev
.
inequality at λ0 = 8.7 × 109 standard deviations.
As mentioned above, an L∞ version of the Stolarsky inequality would eliminate the need for Chebychev inequalities though it might also be very conservative. We do not know whether p(y, ρ̂) has a heavy tailed distribution under
reference distribution 2.
Our permutation points fall into a lattice subset of Rd intersected with the
unit sphere Sd . Our problem of counting the number of such points in a subset
is one that is addressed under the term ‘Geometry of numbers’. According to
a personal communication from Neil Sloane, the standard approach to such
problems is via the volume ratio, which in our setting is p̂1 which does not do
well for our problems.
To have reasonable power to obtain a p-value below  by permutation sampling requires on the order of 1/ permutations each requiring O(n log(n)) computation to generate and O(n) computation to evaluate the inner product. The
cost to compute the standard errors in our method is dominated by a cost
proportional to m3 though there is a very small cost proportional to m4 . In
the range where the first cost dominates, our proposal is advantageous when
m3 = o(n/). Supposing that m0 and m1 are comparable, our advantage holds
when m2 = o(1/). If only the estimate and not the standard error is required,
then our p̂2 and p̂3 cost O(m) once the ρ̂ (cost O(n)) has been computed. Then
the total cost is O(n) compared to the much larger cost O(n/) for sampling.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

26

Acknowledgements
This work was supported by the US National Science Foundation under grants
DMS-1407397 and DMS-1521145. We thank John Robinson and Neil Sloane for
helpful comments.
References
Ackermann, M. and Strimmer, K. (2009). A general modular framework for
gene set enrichment analysis. BMC Bioinformatics, 10:1–20.
Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American Mathematical society, 68(3):337–404.
Barnard, G. A. (1963). Discussion of the spectral analysis of point processes
(by m. s. bartlett). Journal of the Royal Statistical Society, series B, 25:294.
Bilyk, D., Dai, F., and Matzke, R. (2016). Stolarsky principle and energy optimization on the sphere. Technical report.
Blanchet, J. and Glynn, P. (2008). Efficient rare-event simulation for the maximum of heavy-tailed random walks. The Annals of Applied Probability, pages
1351–1378.
Brauchart, J. and Dick, J. (2013). A simple proof of Stolarsky’s invariance
principle. Proceedings of the American Mathematical Society, 141(6):2085–
2096.
Fadista, J., Manning, A. K., Florez, J. C., and Groop, L. (2016). The (in) famous
GWAS p-value threshold revisited and updated for low-frequency variants.
European Journal of Human Genetics, 24(8):1202–1205.
He, H. (2016). Efficient Permutation-Based P-Value Estimation for Gene Set
Tests. PhD thesis, Stanford University.
Jiang, Z. and Gentleman, R. (2007). Extensions to gene set enrichment. Bioinformatics, 23(3):306–313.
Knijnenburg, T. A., Wessels, L. F. A., Reinders, M. J. T., and Shmulevich,
I. (2009). Fewer permutations, more accurate p-values. Bioinformatics,
25(12):i161–i168.
Larson, J. L. and Owen, A. B. (2015). Moment based gene set tests. BMC
Bioinformatics, 16(1):132.
Lee, Y. and Kim, W. C. (2014). Concise formulas for the surface area of the
intersection of two hyperspherical caps. Technical report, Korea advanced
institute of science and technology.
Moran, L. B., Duke, D. C., Deprez, M., Dexter, D. T., Pearce, R. K. B., and
Graeber, M. B. (2006). Whole genome expression profiling of the medial and
lateral substantia nigra in parkinsons disease. Neurogenetics, 7(1):1–11.
Narcowich, F. J., Sun, X., Ward, J. D., and Wu, Z. (2010). Leveque type
inequalities and discrepancy estimates for minimal energy configurations on
spheres. Journal of Approximation Theory, 162(6):1256–1278.
Niederreiter, H. (1992). Random Number Generation and Quasi-Monte Carlo
Methods. SIAM, Philadelphia, PA.

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

27

Reid, N. (1988). Saddlepoint methods and statistical inference. Statistical Science, pages 213–227.
Robinson, J. (1982). Saddlepoint approximations for permutation tests and
confidence intervals. Journal of the Royal statistical society,Series B, pages
91–101.
Scherzer, C. R., AC, A. C. E., Morse, L. J., Liao, Z., Locascio, J. J., Fefer,
D., Schwarzschild, M. A., Schlossmacher, M. G., Hauser, M. A., Vance, J. M.,
Sudarsky, L. R., Standaert, D. G., Growdon, J. H., Jensen, R. V., and Gullans,
S. R. (2007). Molecular markers of early Parkinson’s disease based on gene
expression in blood. Proc Natl Acad Sci, 104(3):955–60.
Stolarsky, K. B. (1973). Sums of distances between points on a sphere. II.
Proceedings of the American Mathematical Society, 41(2):575–582.
Zhang, Y., James, M., Middleton, F. A., and Davis, R. L. (2005). Transcriptional analysis of multiple brain regions in Parkinson’s disease supports the
involvement of specific protein processing, energy metabolism, and signaling
pathways, and suggests novel disease mechanisms. American J Med Genet B
Neuropsychiatry Genet, 137B(1):5–16.
Zhou, C., Wang, H. J., and Wang, Y. M. (2009). Efficient moments-based
permutation tests. In Advances in neural information processing systems,
pages 2277–2285.

11. Appendix
Here we collect up some of the longer proofs.
11.1. Proof of Theorem 3 (Limiting invariance)
Here we show that taking limits as  goes to zero in the formula of Brauchart
and Dick (2013) proves Theorem 3. We use three lemmas, one for each term in
Theorem 2. We use  → 0 as a shorthand for lim1 →0+ lim2 →0+ .
Lemma 4. Let v be defined as in (3.3). Then for ρ̂ ∈ [−1, 1),
Z

1

lim

Z
σd (C(z; t)) −

v (t)

→0

−1

Z
=
Sd

Sd

|p(z, ρ̂) − p̂(ρ̂)|2 dσd (z).

N −1
2
1 X
1C(z;t) (xk ) dσd (z) dt
N
k=0

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

28

Proof. Substituting v we get
Z 1
Z
1
2 + 1(ρ̂ 6 t 6 ρ̂ + 1 )
(p̂(t) − p(z, t))2 dσd (z) dt

d
1
−1
S
Z
Z
1 ρ̂+1
(p̂(t) − p(z, t))2 dσd (z) dt, as 2 → 0+
→
1 ρ̂
Sd
Z
→
(p̂(ρ̂) − p(z, ρ̂))2 dσd (z), as 1 → 0+ .
Sd

Lemma 5. Let v be as in (3.3) with ρ̂ ∈ [−1, 1), and let Kv be given by (3.2).
Then for any x, x0 ∈ Sd ,
lim Kv (x, x0 ) = σd (C(x; ρ̂) ∩ C(x0 ; ρ̂)).

→0

Proof. The argument is essentially the same as for Lemma 4.
Lemma 6. Let v be as in (3.3) with ρ̂ ∈ [−1, 1), and let Kv be given by (3.2).
Then
Z Z
lim
Kv (x, y) dσd (x) dσd (y) = p̂1 (ρ̂)2 .
(11.1)
→0

Sd

Sd

Proof. For any x, y ∈ Sd , the kernel Kv (x, y) is nonnegative and upper bounded
by a constant. Therefore we can take our limit operations
inside the double inR
tegral over x and y. Now lim→0 Kv (x, y) = Sd 1C(z;ρ̂) (x)1C(z;ρ̂) (y) dσd (z).
Therefore the limit in (11.1) is
Z Z Z
1C(z;ρ̂) (x)1C(z;ρ̂) (y) dσd (z) dσd (y) dσd (x) = p̂1 (ρ̂)2
Sd

Sd

Sd

after changing the order of the integrals.
Theorem 3 Let x0 , x1 , . . . , xN ∈ Sd and ρ̂ ∈ [−1, 1]. Then
Z
Sd

|p(z, ρ̂) − p̂1 (ρ̂)|2 dσd (z) =

N −1 N −1
1 X X
σd (C(xk ; ρ̂) ∩ C(x` ; ρ̂)) − p̂1 (ρ̂)2 .
N2
k=0 l=0

Proof. Theorem 2 gives us an identity and applying Lemmas 4, 5 and 6 to both
sides of it establishes (3.4) for ρ ∈ [−1, 1). For ρ̂ = 1 we get the answer by
replacing v by 2 + (1/1 )11−1 6t61 in the lemmas. Replacing z by y and ρ̂ by
t above gives the version in the main body of the article.
11.2. Proof of Lemma 3 (Double inclusion for Model 2)
Proof. We split the proof Rinto four cases and prove them individually. Recall
that P2 (u1 , u2 , u3 , ρ̃, ρ̂) = Sd−1 1(hy, x1 i > ρ̂)1(hy, x2 i > ρ̂) dσd−1 (y ∗ ) where
p
y = ρ̃xc + 1 − ρ̃2 y ∗ .

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

29

Case 1. x1 = x2 = xc , i.e., r1 = r2 = r3 = 0.
Z
P2 (1, 1, 1, ρ̃, ρ̂) =
1(hy, xc i > ρ̂)1(hy, xc i > ρ̂) dσd−1 (y ∗ ) = 1(ρ̃ > ρ̂).
Sd−1

Case 2. x1 = xc 6= x2 , i.e., r1 = 0, r2 > 0, r3 > 0.
Z
P2 (1, u2 , u2 , ρ̃, ρ̂) =
1(hy, xc i > ρ̂)1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
d−1
S
Z
= 1(ρ̃ > ρ̂)
1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
Sd−1

= 1(ρ̃ > ρ̂)P1 (u2 , ρ̃, ρ̂)
where the last step uses Lemma 2.
Case 3. x1 = x2 =
6 xc , i.e., r1 = r2 > 0 = r3 .
Z
P2 (u2 , u2 , 1, ρ̃, ρ̂) =
1(hy, x1 i > ρ̂)1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
d−1
ZS
=
1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
Sd−1

= P1 (u2 , ρ̃, ρ̂).
Case 4. x1 6= x2 6= xc 6= x1 , i.e., r1 , r2 , r3 > 0. We split this case into
subcases. First we assume u2 = −1, so
Z
P2 (u1 , u2 , u3 , ρ̃, ρ̂) =
1(hy, x1 i > ρ̂)1(hy, −xc i > ρ̂) dσd−1 (y ∗ )
d−1
S
Z
= 1(−ρ̃ > ρ̂)
1(hy, x1 i > ρ̂) dσd−1 (y ∗ )
Sd−1

= 1(−ρ̃ > ρ̂)P1 (u1 , ρ̃, ρ̂).
Similarly if u1 = −1, then
P2 (u1 , u2 , u3 , ρ̃, ρ̂) = 1(−ρ̃ > ρ̂)P1 (u2 , ρ̃, ρ̂).
Finally we assume u1 > −1 and u2q> −1, so now |u1 | < 1 and |u2 | <
1. Recall the projections xj = uj cc + 1 − u2j x∗j for j = 1, 2 and introduce
√
further p
projections of y ∗ and x∗2 onto x∗1 : y ∗ = tx∗1 + 1 − t2 y ∗∗ and x∗2 =
∗∗
d
u∗3 x∗1 + 1 − u∗3 x∗∗
and x∗∗
2 belong to a subset of S that is
2 . The residuals y

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

30

isomorphic to Sd−2 . Now we have
P2 (u1 , u2 , u3 , ρ̃, ρ̂)
Z
=
1(hy, x1 i > ρ̂)1(hy, x2 i > ρ̂) dσd−1 (y ∗ )
Sd−1
1

Z
=

Z
q


p
d−1
ωd−2
1 ρ̃u1 + 1 − ρ̃2 1 − u21 t > ρ̂
(1 − t2 ) 2 −1
−1 ωd−1
Sd−2
q
q


p
p
∗∗
∗∗
× 1 ρ̃u2 + 1 − ρ̃2 1 − u22 (tu∗3 + 1 − t2 1 − u∗2
hy
,
x
i)
>
ρ̂
2
3
× dσd−1 (y ∗∗ ) dt

=

1



d−1
ωd−2
ρ̂ − ρ̃u1
p
(1 − t2 ) 2 −1 1 t > p
1 − ρ̃2 1 − u21
−1 ωd−1
Z
q
q


p
p
∗∗
∗∗
×
1 ρ̃u2 + 1 − ρ̃2 1 − u22 (tu∗3 + 1 − t2 1 − u∗2
3 hy , x2 i) > ρ̂

Z

Sd−2

× dσd−1 (y ∗∗ ) dt


1(ρ̃u > ρ̂)1(ρ̃u2 > ρ̂),

R 1 ω1
d−1
d−2
(1 − t2 ) 2 −1 1(t > ρ1 )1(tu∗3 > ρ2 ) dt,
=
−1 ωd−1
R

d−1
1 ωd−2

 −1
(1 − t2 ) 2 −1 1(t > ρ1 )σd−2 (C(x∗∗
2 , √
ω

ρ̃ = ±1
ρ̃ 6= ±1, u∗3 = ±1
ρ2 −tu∗
3

1−t2

d−1

√

1−u∗2
3

)) dt,

ρ̃ 6= ±1, |u∗3 | < 1

where u∗3 , ρ1 , ρ2 are defined in (4.7). Hence, the result follows.
11.3. Proof of Theorem 6 (Second moment under reference
distribution 2)
Proof. Without loss of generality we relabel the values xk so that c = 0. Any
other choice for c is reflected in the number ρe. The second moment is
E(p(y, ρ̂)2 ) =

N −1 N −1
1 X X
P2 (uk , u` , uk,` , ρ̃, ρ̂)
N2

(11.2)

k=0 l=0

where uk,` is obtained via (2.1) from the swap distance rk,` between points xk
and x` . We will partition the sum in (11.2) into the same four cases as in the
proof of Lemma 3.
Case 1, xk = x` = xc , i.e., rk = r` = rk,` = 0. There is only one pair
of (xk , x` ) for this condition. Hence, we get only one term corresponding to
P2 (1, 1, 1, ρ̃, ρ̂) = 1(ρ̃ > ρ̂).
Case 2, xk = xc 6= x` , i.e., rk = 0, r` = rk,` > 0. Consider all pairs of
(xk , x` ) that satisfy this condition and let K2 denote their total contribution

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

31

to (11.2). Then
K2 = 2

N
−1 Z
X
l=1

=2

1(hy, xc i > ρ̂)1(hy, x` i > ρ̂) dσd−1 (y ∗ )

Sd−1

m   
X
m0
m1
r

r=1

r

P2 (1, u(r), u(r), ρ̃, ρ̂).

Case 3, xk = x` 6= xc , i.e., rk = r` > 0 = rk,` . The contribution from terms
of this form is
K3 =

N
−1 Z
X

1(hy, xk i > ρ̂) dσd−1 (y ∗ ) =

Sd−1

k=1

m   
X
m0
m1
r

r=1

r

P1 (u(r), ρ̃, ρ̂).

Case 4, xk 6= x` 6= xc , i.e., rk , r` , rk,` > 0. The contribution of these cases
to the sum is
K4 =

N
−1 N
−1
X
X

Z
Sd−1

k=1 `=1

=

1(hy, xi i > ρ̂)1(hy, xj i > ρ̂) dσd−1 (y ∗ )

1(` 6= k)

X X

X

c(rk , r` , rk,` )P2 (u1 , u2 , u3 , ρ̃, ρ̂).

rk ∈R r` ∈R rk,` ∈R3 (r)

Then the second moment is (1(ρ̃ > ρ̂) + K2 + K3 + K4 )/N 2 .
11.4. Proof of Theorem 7 (Location weighted invariance)
Proof. We follow the technique in Brauchart and Dick (2013). We begin by
showing that Kv,h,x0 as defined in (5.2) is a reproducing kernel. First, Kv,h,x0
is symmetric: Kv,h,x0 (x, y) = Kv,h,x0 (y, x). Next, choose a0 , . . . , aN −1 ∈ R and
NP
−1
x0 , . . . , xN −1 ∈ Sd . Then
ak a` Kv,h,x0 (xk , x` ) equals
k,`=0

Z

1

−1

Z

1

Z

ak a` v(t)h(hz, x0 i)1C(z;t) (xk )1C(z;t) (x` ) dσd (z) dt

Sd k,`=0

Z

=
−1

N
−1
X

Sd

v(t)h(hz, x0 i)

N
−1
X

2

ak 1C(z,t) (xk ) dσd (z) dt

k=0

which is nonnegative. Thus Kv,h,x0 is symmetric and positive definite, and so
by Aronszajn (1950), Kv,h,x0 is a reproducing kernel.
Aronszajn (1950) also shows that a reproducing kernel uniquely defines a
Hilbert space of functions with a specific inner product. Let Hv,h,x0 = H(Kv,h,x0 , Sd )
denote the corresponding reproducing kernel Hilbert space of functions f : Sd →
R with reproducing kernel Kv,h,x0 .

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

32

We now consider functions f1 , f2 : Sd → R which admit the representation
Z 1Z
gi (z; t)1C(z;t) (x) dσd (z) dt, i = 1, 2
(11.3)
fi (x) =
Sd

−1

for functions gi ∈ L2 (Sd ×[−1, 1]). For any fixed y ∈ Sd the function Kv,h,x0 (·, y)
has representation (11.3) via g(z; t) = v(t)h(hz, x0 i)1C(z,t) (y).
For functions with representation (11.3), we define an inner product by
Z
hf1 , f2 iKv,h,x0 =

1

−1

1
v(t)

Z
Sd

1
g1 (z, t)g2 (z, t) dσd (z) dt.
h(hz, x0 i)

(11.4)

For y ∈ Sd and f1 ∈ Hv,h,x0 ,
Z

hf1 , Kv,h,x0 (·, y)iK

v,h,x0

1

Z
1
g1 (z; t)v(t)
=
h(hz, x0 i)1C(z;t) (y) dσd (z) dt
0
−1 v(t) Sd h(hz, x i)
Z 1Z
=
g1 (z, t)1C(z,t) (y) dσd (z) dt
−1

Sd

= f1 (y),
showing that the inner product (11.4) has the reproducing property. By Aronszajn (1950), the inner product in Hv,h,x0 is unique. Functions fi satisfying
(11.3) with hfi , fi iKv,h,x0 < ∞ are in Hv,h,x0 , and (11.4) is the unique inner
product of Hv,h,x0 .
We prove the theorem by equating two different forms of kR(Hv,h,x0 ; · )kKv,h,x0
where
Z
R(Hv,h,x0 ; · ) =

Kv,h,x0 (· , y) dσd (y) −
Sd

N −1
1 X
Kv,h,x0 (· , xk ).
N
k=0

Although R(Hv,h,x0 ; · ) depends on our specific points xi we omit that from the
notation. The reproducing property of Kv,h,x0 yields
hKv,h,x0 (· , xk ), Kv,h,x0 (· , x` )iK

v,h,x0

from which it follows that
Z
Z
0
Kv,h,x (· , y) dσd (y),
Sd

Z

K

v,h,x0


(· , y ) dσd (y )
0

Sd

Z

=

0

Sd

0

Kv,h,x0
0

Kv,h,x0 (y, y ) dσd (y) dσd (y ).
Sd

= Kv,h,x0 (xk , x` )

(11.5)

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

33

Using (11.5) and the linearity of the inner product, we have
hR(Hv,h,x0 ; · ), R(Hv,h,x0 ; · )iK
v,h,x0
Z Z
N −1 Z
2 X
0
0
=
Kv,h,x0 (y, y ) dσd (y) dσd (y ) −
Kv,h,x0 (y, xk ) dσd (y)
N
Sd Sd
Sd
k=0

+

1
N2

N
−1
X

Kv,h,x0 (xk , x` ).

k,`=0

(11.6)
For our second form of kR(Hv,h,x0 ; · )kKv,h,x0 , we write
R(Hv,h,x0 ; · )
Z
N −1
1 X
0
Kv,h,x0 (· , xk )
=
Kv,h,x (· , y) dσd (y) −
N
Sd
k=0
Z

Z 1
Z
N −1
1 X
=
v(t)
1C(z;t) (x)h(hz, xi)
1C(z,t) (y) dσd (y) dt −
1C(z,t) (xk ) dσd (z) dt
N
−1
Sd
Sd
k=0


Z 1
Z
N
−1
1 X
=
v(t)
1C(z,t) (x)h(hz, xi) σd (C(z, t)) −
1C(z,t) (xk ) dσd (z) dt.
N
−1
Sd
k=0

Hence using the definition of the inner product h·, ·iKv,h,x0 , we have
hR(Hv,h,x0 ; x), R(Hv,h,x0 ; x)iK
v,h,x0
Z 1
Z
N −1
2
(11.7)
1 X
1C(x;t) (xk ) dσd (x) dt.
=
v(t)
h(hz, xi) σd (C(x, t)) −
N
−1
Sd
k=0

Combining equations (11.6) and (11.7), we have the generalized location-weighted
version of the Stolarsky invariance principle.
11.5. Proof of Theorem 8 (Spatially weighed invariance)
As in Section 11.1, lim→0 means lim1 →0+ lim2 →0+ and similarly limη→0 denotes limη1 →0+ limη2 →0+ . We prove a series of lemmas first.
Lemma 7. For v (·) and hη (·) defined by equations (3.3) and (5.3),
N −1
2
1 X
1C(z;t) (xk ) dσd (z) dt
lim lim
v (t)
hη (hz, xc i) σd (C(z, t)) −
η→0 →0 −1
N
Sd
k=0
Z
p
∗
2
∗
2
=
|p(ρ̃xc + 1 − ρ̃ y , ρ̂) − p̂1 (ρ̂)| dσd−1 (y ),

Z

1

Z

Sd−1

where p̂1 (ρ̂) = σd (C(y; ρ̂)).

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

34

Proof. This proof is similar to the others. First we take the limit  → 0 yielding
Z
hη (hz, xc i) σd (C(z, ρ̂)) −

lim

η→0

Sd

N
2
1 X
1C(z;ρ̂) (xk ) dσd (z).
N
k=1

√

Making the projection z = sxc + 1 − s2 z ∗ gives
Z 1Z
ωd−1
lim
(1 − s2 )d/2−1 hη (s) ×
η→0 −1 Sd−1 ωd
σd (C(sxc +

p

1 − s2 z ∗ , ρ̂)) −

N
2
1 X
1C(sxc +√1−s2 z∗ ;ρ̂) (xk ) dσd−1 (z ∗ ) ds
N
k=1

Z

p
|p(ρ̂xc + 1 − ρ̂2 y ∗ , ρ̂) − p̂1 (ρ̂)|2 dσd−1 (y ∗ ).

=
Sd−1

Lemma 8. For v (·) and hη (·) defined by equations (3.3) and (5.3) ,
N −1
1 X
Kv ,hη ,xc (xk , x` )
η→0 →0 N 2

lim lim

k,`=0

=

1
N2

N
−1
X

Z

k,`=0

1(hy, xk i > ρ̂)1(hy, x` i > ρ̂) dσd−1 (y ∗ ).

Sd−1

Proof. First, limη→0 lim→0 N −2

NP
−1

Kv ,hη ,xc (xk , x` ) equals

k,`=0

Z
N −1
1 X
lim
hη (hz, xc i)1C(z;ρ̂) (xk )1C(z;ρ̂) (x` ) dσd (z).
η→0 Sd
N2
k,`=0

Projecting z onto xc yields z = sxc +

√

1 − s2 y ∗ and then we have

Z 1
Z
N −1
1 X
ωd−1
2 d/2−1
lim
(1
−
s
)
h
(s)
1C(z;ρ̂) (xk )1C(z;ρ̂) (x` ) dσd−1 (y ∗ )
η
η→0 −1 ωd
N2
d−1
S
k,`=0

N −1 Z
1 X
= 2
1(hy, xk i > ρ̂)1(hy, x` i > ρ̂) dσd−1 (y ∗ ).
N
Sd−1
k,`=0

Lemma 9. For v (·) and hη (·) defined by equations (3.3) and (5.3),
Z Z
lim lim
Kv ,hη ,xc (x, y) dσd (x) dσd (y) = p̂1 (ρ̂)2
η→0 →0

Sd

Sd

Proof. Because Kv ,hη ,xc is nonnegative and uniformly bounded we may take
the limit over  inside the integrals. Now
Z
lim Kv ,hη ,xc (x, y) =
hη (hz, xc i)1C(z;ρ̂) (x)1C(z;ρ̂) (y) dσd (z),
→0

Sd

He, Basu, Zhao & Owen/Stolarsky invariance for permutations

35

and the limit becomes
Z Z Z
lim
hη (hz, xc i)1C(z;ρ̂) (x)1C(z;ρ̂) (y) dσd (z) dσd (x) dσd (y).
η→0

Sd

Sd

Sd

Integrating over z last we get limη→0

R
Sd

hη (hz, xc i)p̂21 (ρ̂) dz = p̂21 (ρ̂).

Lemma 10. Under reference distribution 2
N −1 Z
2 X
lim lim
Kv ,hη ,xc (x, xk ) dσd (x) = 2p̂1 (ρ̂)E(p(y, ρ̂)).
η→0 →0 N
Sd
k=0

Proof. The argument here is similar to the one used for Lemma 9. Take the
limit over  inside the integral and change the order of integration to yield
Z
lim 2p̂1 (ρ̂)

η→0

Sd

N −1
1 X
hη (hz, xc i)1C(z,ρ̂) (xk ) dσd (z).
N
k=0

Substituting the projection z = txc +
Z

1

lim 2p̂1 (ρ̂)

η→0

−1

Z
= 2p̂1 (ρ̂)
Sd−1

√

ωd−1
(1 − t2 )d/2−1 hη (t)
ωd
1
N

N
−1
X

1 − t2 z ∗ produces

Z
Sd−1

N −1
1 X
1C(txc +√1−t2 z∗ ,ρ̂) (xk ) dσd−1 (z ∗ ) dt
N
k=0

1C(ρ̃x +√1−ρ̃2 z∗ ,ρ̂) (xk ) dσd−1 (z ∗ )
c

k=0

= 2p̂1 (ρ̂)E(p(y, ρ̂))
for y under Model 2.
Proof of Theorem 8
Proof. The proof follows from using Lemmas 7 to 10 and Theorem 7.

