Distance Correlation Coefficients for
Lancaster Distributions
arXiv:1502.01413v2 [math.ST] 29 Nov 2016

Johannes Dueck,∗ Dominic Edelmann,† and Donald Richards‡
November 30, 2016

Abstract
We consider the problem of calculating distance correlation coefficients between random vectors whose joint distributions belong to the class of Lancaster
distributions. We derive under mild convergence conditions a general series representation for the distance covariance for these distributions. To illustrate the
general theory, we apply the series representation to derive explicit expressions
for the distance covariance and distance correlation coefficients for the bivariate
normal distribution and its generalizations of Lancaster type, the multivariate
normal distributions, and the bivariate gamma, Poisson, and negative binomial
distributions which are of Lancaster type.
Key words and phrases. Affine invariance; bivariate gamma distribution; bivariate
normal distribution; bivariate negative binomial distribution; bivariate Poisson
distribution; characteristic function; distance correlation coefficient; Lancaster
distributions; multivariate normal distribution.
2010 Mathematics Subject Classification. Primary: 60E05, 62H20; Secondary:
33C05, 42C05, 60E10.
Running head: Distance Correlation and Lancaster Distributions.

1

Introduction

The concepts of distance covariance and distance correlation, introduced by Székely, et
al. [27, 31], have been shown to be widely applicable for measuring dependence between
∗

Institute of Applied Mathematics, University of Heidelberg, Im Neuenheimer Feld 294, 69120
Heidelberg, Germany.
†
Division of Biostatistics, German Cancer Research Center (DKFZ), Im Neuenheimer Feld 280,
69120 Heidelberg, Germany.
‡
Department of Statistics, Pennsylvania State University, University Park, PA 16802, U.S.A.
‡
Corresponding author; e-mail address: richards@stat.psu.edu

1

2

Dueck, Edelmann, and Richards

collections of random variables. As examples of the ubiquity of distance correlation
methods, we note the results on distance correlation given recently by: Székely, et al.
[21, 28, 29, 30, 31], on statistical inference; Sejdinovic, et al. [26], on machine learning;
Kong, et al. [10], on familial relationships and mortality; Zhou [33], on nonlinear time
series; Lyons [17], on abstract metric spaces; Martı́nez-Gomez, et al. [18] and Richards,
et al. [20], on large astrophysical databases; Dueck, et al. [5], on high-dimensional
inference and the analysis of wind data; and Dueck, et al. [6], on a connection with
singular integrals on Euclidean spaces.
A result which is of fundamental importance in distance correlation theory is the
explicit formula for the empirical distance correlation coefficient [31, pp. 2773-2774].
By combining that explicit formula with the fast algorithm of Huo and Székely [9], it
becomes straightforward to apply distance correlation methods to real-world data sets.
On the other hand, the calculation of population distance correlation coefficients
remains an intractable problem generally. Székely, et al. [31, pp. 2785-2786] calculated
the distance correlation coefficient for the bivariate normal distribution; Dueck, et al.
[4, Appendix] extended that result to the general multivariate normal distribution;
and Dueck, et al. [5] calculated the affinely invariant distance correlation coefficient
for the multivariate normal distribution. Otherwise, no such results are yet available
for any other distribution. Hence, the state of distance correlation theory hitherto
is that the empirical coefficients can be calculated readily but the opposite holds for
their population counterparts, generally. Consequently, it was not possible to calculate
distance correlation coefficients explicitly for given nonnormal distributions in terms of
the usual parameters that parametrize these distributions, or to ascertain for nonnormal
distributions any analogs of the limit theorems derived by Dueck, et al. [5, Section 4].
We describe in detail the difficulties arising in attempts to calculate the population
distance correlation coefficients. Let p and q be positive integers. For column vectors
s ∈ Rp and t ∈ Rq , denote by ksk and ktk the standard Euclidean norms on the
corresponding spaces; thus, if s = (s1 , . . . , sp )> then ksk = (s21 + · · · + s2p )1/2 , and
similarly for ktk. Given vectors u and v of the same dimension, we let hu, vi be the
standard Euclidean scalar product of u and v. For jointly distributed random vectors
(X, Y ) ∈ Rp × Rq and non-random vectors (s, t) ∈ Rp × Rq , let


ψX,Y (s, t) = E exp ihs, Xi + iht, Y i ,
√
i = −1, be the joint characteristic function of (X, Y ), and let ψX (s) = ψX,Y (s, 0)
and ψY (t) = ψX,Y (0, t) be the corresponding marginal characteristic functions. For any
z ∈ C, let |z|2 denote the squared modulus of z; also, we use the notation
γp =

π (p+1)/2
.
Γ (p + 1)/2

(1.1)

In the case of distributions with finite first moments, Székely, et al. [31, p. 2772] defined

3

Distance Correlation and Lancaster Distributions
V(X, Y ), the distance covariance between X and Y , to be the positive square-root of
1
V (X, Y ) =
γp γq
2

Z
Rp+q

|ψX,Y (s, t) − ψX (s)ψY (t)|2
ds dt
kskp+1 ktkq+1

(1.2)

and they defined the distance correlation coefficient between X and Y as
V(X, Y )
R(X, Y ) = p
V(X, X)V(Y, Y )

(1.3)

if both V(X, X) and V(Y, Y ) are strictly positive, and otherwise to be zero [31, p. 2773].
For distributions with finite first moments we have 0 ≤ R(X, Y ) ≤ 1, and R(X, Y ) = 0
if and only if X and Y are mutually independent.
For given random vectors X and Y , the fundamental obstacle in calculating the
population distance correlation coefficient (1.3) is the computation of the singular integral (1.2). In particular, the singular nature of the integrand precludes evaluation of
the integral by expanding the numerator, |ψX,Y (s, t) − ψX (s)ψY (t)|2 , and subsequent
term-by-term integration of each of the resulting three terms.
In this paper, we calculate the distance correlation coefficients for pairs (X, Y ) of
random vectors whose joint distributions are in the class of Lancaster distributions, a
class of probability distributions made prominent by Lancaster [15, 16] and Sarmanov
[23]. The distribution functions of the Lancaster family are well-known to have attractive expansions in terms of certain orthogonal functions (Koudou [14]; Diaconis, et
al. [3]). By applying those expansions, we obtain explicit expressions for the distance
covariance and distance correlation coefficients.
Consequently, we derive under mild convergence conditions a general formula for
the distance covariance for the Lancaster distributions. We apply the general formula
to obtain explicit expressions for the distance covariance and distance correlation for
the bivariate normal distributions and some of its generalizations, for the multivariate
normal distributions, and for bivariate gamma, Poisson, and negative binomial distributions. We remark that explicit results can also be obtained for other Lancaster-type
expansions obtained by Bar-Lev, et al. [2]; however, we will omit the details for other
cases because the formulas derived here are entirely representative of other cases.

2

The Lancaster distributions

To recapitulate the class of Lancaster distributions we generally follow the standard
notation in that area, as given by Koudou [13, 14]; cf., Lancaster [16], Pommeret [19],
or Diaconis, et al. [3, Section 6].
Let (X , µ) and (Y, ν) be locally compact, separable probability spaces, such that
2
L (µ) and L2 (ν) are separable. Let σ, a probability measure on X × Y, have marginal

4

Dueck, Edelmann, and Richards

distributions µ and ν; then there exist functions Kσ and Lσ such that
σ(dx, dy) = Kσ (x, dy)µ(dx) = Lσ (dx, y)ν(dy).
We note that Kσ and Lσ represent the conditional distributions of Y given X = x, and
X given Y = y, respectively.
Let C denote a countable index set with a zero element, denoted by 0. Let {Pn : n ∈
C} and {Qn : n ∈ C} be sequences of functions on X and Y which form orthonormal
bases for the separable Hilbert spaces L2 (µ) and L2 (ν), respectively. We assume, by
convention, that P0 ≡ 1 and Q0 ≡ 1.
Since the tensor product Hilbert space L2 (µ ⊗ ν) ≡ L2 (µ) ⊗ L2 (ν) is separable there
holds, for σ ∈ L2 (µ ⊗ ν), the expansion
XX
σ(dx, dy) =
ρm,n Pm (x)Qn (y) µ(dx) ν(dy),
(2.1)
m∈C n∈C

(x, y) ∈ X × Y. Letting δm,n denote Kronecker’s delta, the probability measure σ is
called a Lancaster distribution if there exists a nonnegative sequence {ρn : n ∈ C} such
that
Z
Pm (x) Qn (y) σ(dx, dy) = ρm δm,n
for all m, n ∈ C; in particular, ρ0 = 1. The sequence {ρn : n ∈ C} is called a Lancaster
sequence, and the expansion (2.1) reduces to
X
ρn Pn (x)Qn (y)µ(dx)ν(dy).
σ(dx, dy) =
n∈C

Koudou [13, pp. 255–256] characterized the Lancaster sequences {ρn : n ∈ C} such that
the associated probability distribution σ is absolutely continuous with respect to µ ⊗ ν
and has Radon-Nikodym derivative
X
σ(dx, dy)
=
ρn Pn (x) Qn (y) ∈ L2 (µ ⊗ ν),
µ(dx) ν(dy) n∈C
(x, y) ∈ X × Y.
In the sequel, we consider the case in which X = Rp and Y = Rq and the underlying
random vectors X ∈ Rp and Y ∈ Rq have joint distribution σ and marginal distributions
µ and ν, respectively. We assume that µ, ν, and σ are absolutely continuous with
respect to Lebesgue measure or counting measure on the respective sample spaces
and we denote their corresponding probability density functions by φX , φY , and φX,Y ,
respectively. This yields the expansion,
X
φX,Y (x, y) = φX (x) φY (y)
ρn Pn (x) Qn (y).
(2.2)
n∈C

We will refer to (2.2) as the Lancaster expansion of the joint density function φX,Y .

5

Distance Correlation and Lancaster Distributions

3

Examples of Lancaster expansions

In this section, we provide examples of Lancaster expansions (2.2) for the bivariate normal distribution and some of its generalizations, the multivariate normal distributions,
and the bivariate gamma, Poisson, and negative binomial distributions. In the sequel,
we denote by N0 the set of nonnegative integers.

3.1

The bivariate normal distribution and some of its generalizations

Let (X, Y ) follow a bivariate normal distribution with mean vector 0 and covariance
matrix


1 ρ
Σ=
,
ρ 1
denoted by (X, Y ) ∼ N2 (0, Σ). The joint probability density function of (X, Y ) is
 2

1
1
x + y 2 − 2ρ x y
2 −2
φX,Y (x, y; ρ) =
(1 − ρ ) exp −
,
2π
2(1 − ρ2 )
x, y ∈ R, and the marginal density functions are given by

1
φX (x) = φY (x) = √ exp − 21 x2 .
2π
In this case, the index set C is N0 . For n ∈ N0 , let
 d n

exp − 12 x2 ,
Hn (x) = (−1)n exp 21 x2 )
dx
x ∈ R, denote the nth Hermite polynomial, n = 0, 1, 2, . . .. It is well-known that
the polynomials {Hn : n ∈ N0 } are orthogonal with respect to the standard normal
distribution and form a complete orthogonal basis for the Hilbert space L2 (X). Also,
the Lancaster expansion of φX,Y is given by the classical formula of Mehler, which states
that, for x, y ∈ R,
φX,Y (x, y; ρ) = φX (x) φY (y)

∞
X
ρn
n=0

n!

Hn (x) Hn (y),

(3.1)

and this series converges absolutely for all x, y ∈ R.
We remark that there are numerous extensions of Mehler’s formula which represent Lancaster-type expansions for generalizations of the bivariate normal distribution.
Sarmanov and Bratoeva [25] consider series expansions of the form
φX,Y (x, y) = φX (x) φY (y)

∞
X
ρn
n=0

n!

Hn (x)Hn (y),

(3.2)

6

Dueck, Edelmann, and Richards

P
2
x, y ∈ R, where the sequence of real numbers {ρn : n = 0, 1, 2, . . .} satisfies ∞
n=0 ρn <
∞. Sarmanov and Bratoeva proved that for the expansion (3.2) to be nonnegative, and
therefore to be a valid probability density function, it is necessary and sufficient that
the sequence {ρn } be the moment sequence of a random variable ξ supported on the
interval [−1, 1].
An example of this generalization is the case in which φX,Y is a mixture of bivariate
normal densities; put in the formula
1
2



φX,Y (x, y; ρ) + φX,Y (x, y; −ρ)
∞
X
ρ2n
= φX (x)φY (y)
H2n (x)H2n (y);
(2n)!
n=0

φX,Y (x, y) =

(3.3)

this corresponds to the case in (3.2) in which
ρn =

(
ρn , n even
0,

n odd

This mixture density also provides an example of a distribution for which the Pearson
correlation coefficent is zero whereas the distance correlation is positive.

3.2

The multivariate normal distribution

Let X ∈ Rp and Y ∈ Rq be random vectors such that (X, Y ) ∼ Np+q (0, Σ), a (p + q)dimensional multivariate normal distribution with mean vector 0 and positive definite
covariance matrix


ΣX ΣXY
Σ=
(3.4)
ΣYX ΣY
where ΣX , ΣY , and ΣXY = ΣYX 0 are p × p, q × q and p × q matrices, respectively. We
denote by φX,Y the joint probability density function of (X, Y ), and by φX and φY the
marginal density functions of X and Y , respectively.
We now describe the Lancaster expansion of φX,Y , a result derived in [32]. In this
case, the index set C is Np×q
0 , the set of p × q matrices with nonnegative integer entries.
Q Q
For a matrix of summation indices N = (Nrc ) ∈ N0p×q , define N ! = pr=1 qc=1 Nrc !.
For r = 1, . . . , p, let
q
X
Nr q =
Nrc
c=1

and set N ∗ q = (N 1 q , . . . , N p q ). Similarly, for each c = 1, . . . , q, define
N qc =

p
X
r=1

Nrc

7

Distance Correlation and Lancaster Distributions
and set N q∗ = (N q1 , . . . , N qq ). Further, we define
N qq =

p
q
X
X

Nrc ,

r=1 c=1

P
P
and note that N q q = pr=1 N r q = qc=1 N qc .
Denoting by (ΣXY )rc the (r, c)th entry of ΣXY , we also define
ΣXY

N

p
q
Y
Y
[(ΣXY )rc ]Nrc .
=
r=1 c=1

We now introduce the multivariate Hermite polynomials. For any p ∈ N, k =
k
(k1 , . . . , kp ) ∈ Np0 , and x = (x1 , . . . , xp ) ∈ Rp , define xk = xk11 · · · xpp and define the
differential operator,

k 
k1

kp
∂
∂
∂
−
= −
··· −
.
∂x
∂x1
∂xp
The kth multivariate Hermite polynomial with respect to the marginal density function
φX is defined as

k
1
∂
Hk (x; ΣX ) =
−
φX (x).
(3.5)
φX (x)
∂x
The Lancaster expansion of the multivariate normal density function φX,Y is given by
the generalized Mehler formula [32]:
φX,Y (x, y) = φX (x) φY (y)

X ΣXY N
HN ∗ q (x; ΣX ) HN q∗ (y; ΣY ),
N!
p×q

(3.6)

N ∈N0

with absolute convergence for all x ∈ Rp , y ∈ Rq .
To calculate the affinely invariant distance correlation coefficient between X and Y ,
as defined by Dueck, et al. (2014), we need the Lancaster expansion of the joint density
e = ΣX −1/2 X and Ye = ΣY −1/2 Y . It is
function of the standardized random vectors X
e Ye ) ∼ Np+q (0, Λ) where
straightforward to verify that (X,

Λ=

Ip
ΛXY
0
ΛXY
Iq


(3.7)

with ΛXY = ΣX −1/2 ΣXY ΣY −1/2 , and then we deduce from (3.6) that the Lancaster
e Ye ) is
expansion for (X,
φX,
e Ye (x, y) = φX
e (x) φYe (y)

X ΛXY N
HN ∗ q (x; Ip ) HN q∗ (y; Iq ).
N!
p×q

N ∈N0

(3.8)

8

3.3

Dueck, Edelmann, and Richards

The bivariate gamma distribution

The Lancaster expansion for a bivariate gamma distribution, which was derived by
Sarmanov [24, 22], can be stated as follows (see Kotz, et al. [11, pp. 437–438]).
For α > −1 and n ∈ N0 , the classical Laguerre polynomial is defined by
 d n
1 −α
(α)
e
Ln (x) = x exp(x)
xn+α exp(−x)
n!
dx
n
(α + 1)n X (−n)j xj
=
,
n!
(α + 1)j j!
j=0

(3.9)

x > 0, where
(α)n =

Γ(α + n)
= α(α + 1) · · · (α + n − 1),
Γ(α)

n = 0, 1, 2, . . ., denotes the rising factorial. By standardizing the classical Laguerre
polynomial, we obtain the orthonormal version [8],
L(α)
n (x)


=

(α + 1)n
n!

−1/2

e(α) (x)
L
n


=

(α + 1)n
n!

1/2 X
n
(−n)j xj
.
(α + 1)j j!
j=0

Let λ ∈ (0, 1), and let α and β satisfy α ≥ β > 0. Sarmanov [24, 22] derived for
certain bivariate gamma random variables (X, Y ) the joint probability density function,
φX,Y (x, y) = φX (x) φY (y)

∞
X

an L(α−1)
(x) L(β−1)
(y),
n
n

(3.10)

n=0

x, y > 0, where


(β)n
an =
(α)n

1/2

λn ,

(3.11)

n = 0, 1, 2, . . .. The corresponding marginal density functions are
φX (x) =

1
xα−1 exp(−x)
Γ(α)

φY (y) =

1
y β−1 exp(−y),
Γ(β)

and

which we recognize as the density functions of one-dimensional gamma random variables
with index parameters α and β, respectively.
We remark that if α = β then the density function (3.10) reduces to the KibbleMoran bivariate gamma density function, Corr(X, Y ) = λ [11, pp. 436–437], and
(3.10) represents the Lancaster expansion for (X, Y ). On the other hand, if α 6= β then
Corr(X, Y ) 6≡ λ.

9

Distance Correlation and Lancaster Distributions
More generally, Griffiths [8] showed that a series expansion of the form
φX,Y (x, y) = φX (x) φY (y)

∞
X

(y)
(x) L(β−1)
ρn L(α−1)
n
n

(3.12)

n=0

represents a valid bivariate probability density if and only if

1/2
(β)n
ρn =
λn ,
(α)n

(3.13)

where λn is the moment sequence of a random variable ξ concentrated on [0, 1].

3.4

The bivariate Poisson distribution

For a > 0 and x, n ∈ N0 , let
Cn (x; a) =

n
 an 1/2 X

n!

  
n x k!
(−1)
k
k ak
k=0
k

(3.14)

denote the Poisson-Charlier polynomial of degree n. For λ ∈ [0, 1], Koudou [14, Section
5] (cf., Bar-Lev, et al. [2], Pommeret [19]) showed that there exists a bivariate random
vector (X, Y ) with probability density function
φX,Y (x, y) = φX (x) φY (y)

∞
X

λn Cn (x; a) Cn (y; a),

(3.15)

n=0

x, y ∈ N0 . The corresponding marginal density functions φX and φY are given by
ak exp(−a)
,
k!
k ∈ N0 , so that X and Y are distributed marginally according to a Poisson distribution
with parameter a. The series (3.15) is an expansion of Lancaster type, a special case
of (2.2), and the resulting distribution is called a bivariate Poisson distribution.
φX (k) = φY (k) =

3.5

The bivariate negative binomial distribution

The orthonormal polynomials for the classical univariate negative binomial distribution
are the (normalized) Meixner polynomials, given by
n
 cn (β) 1/2 X
k
(−n)k (−x)k
n
β,c
1 − c−1 ,
(3.16)
Mn (x) =
n!
(β)
k k!
k=0
for β > 0, 0 < c < 1, and x ∈ N. Koudou [14, Section 6] showed, by an approach
similar to that used for the bivariate Poisson distribution, that there exists a bivariate
random variable (X, Y ), with identical marginal negative binomial densities,
φX (x) = φY (x) = (1 − c)β

cx (β)x
,
x!

10

Dueck, Edelmann, and Richards

x ∈ N0 , and with joint probability density function,
φX,Y (x, y) = φX (x) φY (y)

∞
X

λn Mnβ,c (x) Mnβ,c (y),

(3.17)

n=0

where x, y ∈ N0 , and 0 ≤ λ < 1. The expansion (3.17) represents a Lancaster expansion
of the joint density function.

4

Distance correlation coefficients for Lancaster distributions

In this section, we derive a general series expression for the distance correlation coefficients for Lancaster distributions with density functions of the form (2.2). For a joint
density function φX,Y given by (2.2) and n ∈ C, we introduce the notation
Pn (s) = E exp(i hs, Xi) Pn (X),

(4.1)

Qn (t) = E exp(i ht, Y i) Qn (Y ),

(4.2)

s ∈ Rp , and
t ∈ Rq . To verify that the expectation in (4.1) converges absolutely for all s ∈ Rp , we
apply the Cauchy-Schwarz inequality to obtain
1/2
1/2
E| exp(ihs, Xi)Pn (X)| ≤ E| exp(ihs, Xi)|2
· E|Pn (X)|2
= 1,
because {Pn : n ∈ C} is an orthonormal basis for the Hilbert space L2 (µ). In particular,
|Pn (s)| ≤ E| exp(ihs, Xi)Pn (X)| ≤ 1,
for all s ∈ Rp and, similarly, |Qn (t)| ≤ 1 for all t ∈ Rq .
In the following result, we will use the notation
Z
ds
Aj,k =
Pj (s) Pk (−s)
kskp+1
Rp
and

Z

dt
,
ktkq+1
Rq
j, k ∈ C, whenever these integrals converge absolutely.
We now state the main result.
Bj,k =

Qj (t) Qk (−t)

Theorem 4.1. Suppose that the random vectors X ∈ Rp and Y ∈ Rq have the joint
probability density function (2.2). Then,
X
X
1
V 2 (X, Y ) =
ρj ρk Aj,k Bj,k ,
(4.3)
γp γq j∈C,j6=0 k∈C,k6=0
whenever the sum converges absolutely.

Distance Correlation and Lancaster Distributions

11

Proof. Rewriting the Lancaster expansion (2.2) in the form,
X
φX,Y (x, y) − φX (x) φY (y) = φX (x) φY (y)
ρn Pn (x) Qn (y),
n∈C,n6=0

and taking Fourier transforms on both sides of this identity, we obtain for all s ∈ Rp
and t ∈ Rq the expansion
X
ψX,Y (s, t) − ψX (s) ψY (t) =
ρn Pn (s) Qn (t).
(4.4)
n∈C,n6=0

This identity is valid subject to the requirement that we may interchange summation
and integration, which is justified by the assumption that the sum in the final result
converges absolutely. Using (4.4) we deduce that


|ψX,Y (s, t) − ψX (s)ψY (t)|2 = ψX,Y (s, t) − ψX (s)ψY (t) ψX,Y (s, t) − ψX (s)ψY (t)
X
X
=
ρj ρk Pj (s) Pk (−s) Qj (t) Qk (−t).
j∈C,j6=0 k∈C,k6=0

Next, we integrate this expansion with respect to the measures ds/kskp+1 and dt/ktkq+1 ;
this requires that we again interchange summation and integration which, by assumption, we are able to do. On carrying through these procedures, we obtain (4.3).
Remark 4.2. We note that Theorem 4.1 can be extended to the more general αdistance covariance and correlation measures treated by Székely, et al. [31, p. 2784].
For α ∈ (0, 2), define

2 π p/2 Γ 1 − 21 α

γp,α =
α2α Γ 21 (p + α)
and let

Z
1
|ψX,Y (s, t) − ψX (s)ψY (t)|2
)=
ds dt
γp,α γq,α Rp+q
kskp+α ktkq+α
be the α-distance covariance between the random vectors X and Y . Further, define
Z
ds
Aj,k (α) =
Pj (s) Pk (−s)
kskp+α
Rp
Vα2 (X, Y

and

Z

dt
,
ktkq+α
Rq
j, k ∈ C, whenever these integrals converge absolutely. Then, the extension of Theorem
4.1 to the α-distance covariance measures is that
X
X
1
Vα2 (X, Y ) =
ρj ρk Aj,k (α) Bj,k (α),
γp,α γq,α j∈C,j6=0 k∈C,k6=0
Bj,k (α) =

Qj (t) Qk (−t)

whenever the sum converges absolutely. The proof of this result is similar to the proof
of Theorem 4.1.

12

5

Dueck, Edelmann, and Richards

Examples

In this section, we establish the versatility of Theorem 4.1 by applying it to compute
the distance correlation coefficients for the bivariate normal, multivariate normal, and
bivariate gamma, Poisson, and negative binomial distributions. We verify for each
example the absolute convergence of the series resulting from Theorem 4.1, for that
convergence property cannot be obtained in general from the general theorem. In
developing each example, we retain the corresponding notation in Section 3. We also
remark that the crucial singular integral in the theory of distance correlation [6, 31] is
evaluated in terms of the gamma function; however, even slight generalizations of that
integral can be evaluated only in terms of the Gaussian or the confluent hypergeometric
series; this explains the appearance of those series in the ensuing examples.

5.1

The bivariate normal distribution and some of its generalizations

In the sequel, we use the standard double-factorial notation,


if n = −1, 0

1,
n!! = n(n − 2)(n − 4) · · · =

n(n − 2)(n − 4) · · · 1,


n(n − 2)(n − 4) · · · 2,

if n = 1, 3, 5, 7, . . .

if n = 2, 4, 6, 8, . . .

Proposition 5.1. Let (X, Y ) ∼ N2 (0, Σ), a bivariate normal distribution with correlation coefficient ρ. Then,
2

V (X, Y ) = 4π

−1

∞
X
((2` − 3)!!)2
`=1

(2`)!


1 − 2−(2`−1) ρ2` ,

(5.1)

and this series converges absolutely for all ρ ∈ (−1, 1).
Proof. Starting with the Lancaster expansion of the bivariate normal density function, as given in (3.1), and using the definitions of Pn and Qn in (4.1) and (4.2),
respectively, we obtain by substitution and integration-by-parts,
Z ∞

1
Pn (s) = Qn (s) =
exp(isx) √ exp − 12 x2 Hn (x)dx
2π
−∞

= (is)n exp − 21 s2 ,
s ∈ R. Therefore,
k j+k

Z

∞

Aj,k = Bj,k = (−1) i
sj+k−2 exp(−s2 ) ds,
−∞
(

k j+k 1/2 1 (j+k−2)/2
(−1) i
π
(j + k − 3)!!, if j + k is even
2
=
0,
otherwise

13

Distance Correlation and Lancaster Distributions
since the latter integral is a moment of the N (0, 21 ) distribution. By Theorem 4.1, we
obtain
2
4 X ρj+k 1 j+k
V 2 (X, Y ) =
(j
+
k
−
3)!!
.
π j, k > 0 j! k! 2
j+k even

Setting j + k = 2` with ` ≥ 1, the double series reduces to
∞
X 1
4 X 2` 1 2`
V (X, Y ) =
ρ ( 2 ) ((2` − 3)!!)2
π `=1
j! k!
j,k≥1
2

j+k=2`

4
=
π
=

4
π

∞
X
`=1
∞
X

2`−1
− 3)!!)2 X
(2`)!
(2`)!
j! (2` − j)!
j=1

((2`
ρ2` ( 12 )2`
ρ2` ( 12 )2`

`=1

((2` − 3)!!)2 2`
(2 − 2),
(2`)!

which is the same as (5.1).
The absolute convergence of (5.1) can be verified by comparison with a geometric
series. Moreover, it can be shown that the series reduces to the explicit formula,

p
p
4
ρ sin−1 ρ + 1 − ρ2 − ρ sin−1 ρ/2 − 4 − ρ2 + 1
π

(5.2)

which is identical with the result obtained by Székely, et al. [31, pp. 2785-2786].
Having obtained V(X, Y ), we let ρ → 1− to obtain the distance variances V(X, X)
and V(Y, Y ); here, we are applying a well-known result that if (X, Y ) ∼ N2 (0, Σ) where
Var(X) = Var(Y ) and ρ = 1 then X = Y , almost surely. By applying properties of
Gauss’ hypergeometric series, as was done by Dueck, et al. [5, p. 2318], we obtain
√
4 4( 3 − 1)
2
2
.
V (X, X) = V (Y, Y ) = −
3
π
It is straightforward to extend the above results to generalizations of the type given
in equation (3.2).
Corollary 5.2. Let (X, Y ) be a bivariate random variable distributed according to a
density function as given in (3.2). Then
V 2 (X, Y ) =

4
π

X
j, k > 0
j+k even

ρj ρk
j! k!


1 j+k
2

2
(j + k − 3)!! .

(5.3)

For the example given in (3.3), the series expansion in (5.3) reduces to an explicit
formula similar to (5.2).

14

Dueck, Edelmann, and Richards

Corollary 5.3. Let (X, Y ) be a bivariate random variable distributed according to a
density function as given in (3.2) with ρn = ρn for n even and ρn = 0 for n odd. Then
p
p
4  ρ −1
3
2
−1
1
2
2
V (X, Y ) =
sin ρ + 2 1 − ρ − ρ sin ρ/2 − 4 − ρ +
.
(5.4)
π 2
2
Proof. Proceeding as in the case of the bivariate normal distribution, we obtain
∞
X
1
4 X 2` 1 2`
V (X, Y ) =
ρ ( 2 ) ((2` − 3)!!)2
π `=1
(2j)! (2k)!
j,k≥1
2

j+k=`

4
=
π

∞
X

`−1
− 3)!!)2 X
(2`)!
(2`)!
(2j)! (2` − 2j)!
j=1

((2`
ρ2` ( 12 )2`

`=1

∞
4 X 2` 1 2` ((2` − 3)!!)2 2`−1
ρ (2)
(2
− 2).
=
π `=1
(2`)!

Using the standard notation, 2 F1 , for Gauss’ hypergeometric function we see that
∞
4 X 2` 1 2` (2` (− 12 )` )2 2`−1
(2
− 2)
V (X, Y ) =
ρ (2)
π `=1
22` `! ( 12 )`

i
4 h1
2
1
1
1 1
1 1 1 2
=
F
(−
F
(−
.
,
−
;
;
ρ
)
−
1
−
2
,
−
;
;
ρ
)
−
1
2 1
2 1
2
2 2
2
2 2 4
π 2
It is well-known (see Andrews, Askey, and Roy [1, , pages 64 and 94]) that
2

2
1
1 1
2 F1 (− 2 , − 2 ; 2 ; ρ )

= ρ sin−1 ρ + (1 − ρ2 )1/2 .

On applying this formula to the above expression, we obtain (5.4).

5.2

The multivariate normal distribution

In this subsection, we will make extensive use of the notation N r q , N qc , N ∗ q , N q∗ ,
and N q q from Subsection 3.2 for the multi-index matrix N ∈ Np×q
0 . We now establish
the following result.
Proposition 5.4. Suppose that (X, Y ) ∼ Np+q (0, Σ), where Σ is given in (3.4). Then
e 2 (X, Y ), is given by
the affinely invariant distance covariance, V
e 2 (X, Y ) =
V

1
γp γq

X
J6=0,K6=0

AJ,K BJ,K

ΛJXY ΛK
XY
,
J ! K!

(5.5)

where the sums are taken over all non-zero J , K ∈ Np×q
such that all components of
0
J ∗ q + K ∗ q and J q∗ + K q∗ are even,
 p

Γ 21 (J q q + K q q − 1) Y 1

AJ,K =
Γ 2 (J r q + K r q + 1)
(5.6)
1
1
Γ 2 (J q q + K q q ) + 2 p r=1

15

Distance Correlation and Lancaster Distributions
and
BJ,K

 q
Γ 12 (J q q + K q q − 1) Y

Γ
=
Γ 12 (J q q + K q q ) + 21 q c=1

1
(J qc
2


+ K qc + 1) .

(5.7)

Proof. In this case, the index set C is Np×q
0 , and we write the Lancaster expansion
e
e
(3.8) of (X, Y ) in the form
φX,
e Ye (x, y) − φX
e (x) φYe (y) = φX
e (x) φYe (y)

X ΛN
XY
HN ∗ q (x; Ip ) HN q∗ (y; Iq ).
N!
N 6=0

e we apply the definition
To calculate the Fourier transform PN corresponding to X,
(3.5) of the multivariate Hermite polynomials and integration-by-parts to deduce that
for s ∈ Rp ,
Z
PN (s) =
exp(ihs, xi) φXe (x) HN ∗ q (x; Ip ) dx
Rp
 N ∗ q
Z
∂
N qq
= (−1)
exp(ihs, xi)
φXe (x) dx
∂x
Rp
 N ∗ q
Z
∂
φXe (x)
=
exp(ihs, xi) dx
∂x
Rp
Z
N∗ q
φXe (x) exp(ihs, xi) dx
= (is)
=i

N

q q N∗ q
s

Rp

exp(− 21 hs, si).

Similarly,
QN (t) = i N q q tN q∗ exp(− 12 ht, ti),
t ∈ Rq . Therefore,
Z
PJ (s) PK (−s)
Rp

ds
= (−1)K q q i J q q +K q q
p+1
ksk

Z

sJ ∗ q +K ∗ q exp(−hs, si)

Rp

ds
.
kskp+1

We now change variables to hyperspherical coordinates: s = rω, where r > 0 and
ω = (ω1 , . . . , ωp ) ∈ S p−1 , the unit sphere in Rp . Then the latter integral reduces to
Z
Z
J q q +K q q −2
2
r
exp(−r ) dr ·
ω J ∗ q +K ∗ q dω.
S p−1

R+

The integral over R+ is evaluated by replacing r by r1/2 , and we obtain its value as

1
1
Γ
(J
+
K
−
1)
.
q
q
q
q
2
2
It is easy to see that the integral over S p−1 equals zero if any component of J ∗ q +K ∗ q
is odd. For the case in which each component of J ∗ q + K ∗ q is even, we obtain
Z
ω J ∗ q +K ∗ q dω = A(S p−1 ) E(ω J ∗ q +K ∗ q ),
S p−1

16

Dueck, Edelmann, and Richards

where A(S p−1 ) = 2π p/2 /Γ( 12 p) is the surface area of S p−1 and ω now is a uniformly distributed random vector on S p−1 . It is well-known that the random vector (ω12 , . . . , ωp2 ) ∼
D( 21 , . . . , 12 ), a Dirichlet distribution with parameters ( 12 , . . . , 12 ); so, by a classical formula for the moments of the Dirichlet distribution [11, p. 488],
Q
Γ( 12 p) pr=1 Γ( 12 (J r q + K r q + 1))
J ∗ q +K ∗ q
)=
E(ω
.
[Γ( 12 )]p Γ( 12 (J q q + K q q ) + 21 p)
Collecting together these results, we obtain
Z
ds
PJ (s) PK (−s)
= (−1)K q q (−1)(J q q +K q q )/2 AJ,K ,
p+1
ksk
Rp
where AJ,K is given in (5.6). A similar expression can be obtained for
Z
dt
QJ (t) QK (−t) q+1 ,
ktk
Rq
from which the final result (5.5) follows.
As a consequence of Proposition 5.4, we now derive the value of the affinely ine 2 (X, X), when X has a multivariate normal distribution.
variant distance variance, V
We remark that the derivation of this result given by Dueck, et al. [5, Corollary 3.3]
utilized the theory of zonal polynomials, whereas the proof which we now give is by a
simpler method.
Corollary 5.5. (Dueck, et al. [5, Corollary 3.3]) Suppose that (X, Y ) ∼ N2p (0, Σ),
where Σ in (3.4) satisfies ΛXY = ρ Ip in (3.7). Then
e 2 (X, Y ) = 4π
V

2


γp−1
2
1 1
1
1 1
1 2
1
,
−
;
p;
ρ
)
−
2
F
(−
,
−
;
p;
ρ
)
+
1
F
(−
2
1
2
1
2
2 2
2
2 2
4
γp2

(5.8)

and
γ2
e 2 (X, X) = 4π p−1
V
γp2

"

#

Γ( 21 p) Γ( 21 p + 1)
 1
2 − 2 2 F1 − 12 , − 12 ; 12 p; 41 + 1 .
Γ 2 (p + 1)

(5.9)

Proof. Setting p = q and ΛXY = ρIp in (5.5), we obtain
e 2 (X, Y ) =
γp2 V

XX
J6=0,K6=0
J ∗ q +K ∗ q ,J q∗ +K

AJ,K BJ,K

(ρIp )J+K
.
J ! K!

q∗ even

By decomposing the set of all (J , K) into a disjoint union,
{(J , K) : J + K 6= 0} = {(J , K) : J 6= 0, K 6= 0}
∪ {(J , K) : J = 0, K 6= 0}
∪ {(J , K) : J 6= 0, K = 0},

(5.10)

17

Distance Correlation and Lancaster Distributions
and noting that the summand in (5.10) is symmetric in (J , K), we obtain
XX

e 2 (X, Y ) =
γp2 V

J+K6=0
J ∗ q +K ∗ q ,J q∗ +K

AJ,K BJ,K

(ρIp )J+K
J ! K!

q∗ even

XX

−2

AJ,K BJ,K

J=0,K6=0
J ∗ q +K ∗ q ,J q∗ +K

(ρIp )J+K
.
J ! K!

(5.11)

q∗ even

Note that
(ρIp )

J+K

p
p
p
Y
Y
Y
Y
Jrc +Krc
(ρδrr )Jrr +Krr ·
(ρδrc )
=
(ρδrc )Jrc +Krc .
=
r=1

r=1 c=1

r6=c

This is non-zero iff Jrc = Krc = 0 for all r 6= c, in which case J and K are diagonal
matrices, and then we have
(ρIp )J+K = ρJ11 +···+Jpp +K11 +···+Kpp ,
J ∗ q + K ∗ q = (J11 + K11 , . . . , Jpp + Kpp ) = J q∗ + K q∗ ,
and
J! =

p
Y

Jrr !,

K! =

r=1

p
Y

Krr !.

r=1

Therefore,
AJ,K = BJ,K

 p
Γ 21 (J11 + · · · + Jpp + K11 + · · · + Kpp − 1) Y

=
Γ
Γ 21 (J11 + · · · + Jpp + K11 + · · · + Kpp ) + 21 p r=1

1
(Jrr
2


+ Krr + 1) ,

and this yields for the first term in (5.11)
XX
J+K6=0
J ∗ q +K ∗ q ,J q∗ +K

=

AJ,K BJ,K

(ρIp )J+K
J ! K!

q∗ even

X

X

A2J,K

n1 ,...,np ≥0 J11 +K11 =2n1 ,...,Jpp +Kpp =2np
n1 +···+np 6=0

=

X

ρ2n1 +···+2np

n1 ,...,np ≥0 and even
n1 +···+np 6=0

×



ρJ11 +···+Jpp +K11 +···+Kpp
J ! K!

 p

 2
Γ n1 + · · · + np − 21 Y
1

Γ nr + 2
Γ n1 + · · · + np + 12 p r=1

X
J11 +K11 =2n1 ,...,Jpp +Kpp =2np

1
.
J ! K!

(5.12)

18

Dueck, Edelmann, and Richards

Since
p

X
J11 +K11 =2n1 ,...,Jpp +Kpp =2np

Y
1
=
J ! K! r=1



X
Jrr +Krr =2nr

 Y
p
1
22nr
=
,
Jrr ! Krr !
(2nr )!
r=1

we obtain that (5.12) equals


X
n1 ,...,np ≥0
n1 +···+np 6=0

 p
 p 2nr 2nr
 2Y
Γ n1 + · · · + np − 12 Y
2 ρ
1

Γ nr + 2
1
(2nr )!
Γ n1 + · · · + np + 2 p r=1
r=1


p 
∞ 
1 2
X
X
Y
Γ n − 21 2
Γ
n
+
r
2
 (2ρ)2n
=
.
1
(2n
)!
Γ
n
+
p
r
2
n=1
n +···+np =n r=1
1

However,
X
n1 +···+np

2
p 
Y
X
Γ nr + 12
=
(2nr )!
=n r=1
n +···+n

p

1

p

π
= 2n
2

2
p  1
Y
( 2 )nr Γ( 12 )
22nr nr !( 12 )nr
=n r=1

X

p
Y
( 1 )nr

n1 +···+np =n r=1

2

nr !

π ( 12 p)n
,
22n n!
p

=

so we obtain that the first term in (5.11) is given by


∞ 
∞ 
p 1
X
X
Γ n − 12 2
Γ(− 12 )(− 12 )n 2 2n p ( 21 p)n
2n π ( 2 p)n
 (2ρ) 2n
=
ρ π
1
2
n!
n!
Γ( 12 p)( 12 p)n
Γ
n
+
p
2
n=1
n=1
∞
[Γ(− 12 )]2 X (− 12 )n (− 21 )n ρ2n
n!
[Γ( 12 p)]2 n=1
( 12 p)n


2
2
1
1 1
= 4πγp−1
2 F1 (− 2 , − 2 ; 2 p; ρ ) − 1 .

= πp

By similar arguments, we obtain for the second term in (5.11):
XX
J=0,K6=0
J ∗ q +K ∗ q ,J q∗ +K q∗ are even

AJ,K BJ,K



(ρIp )J+K
2
1
1 1
1 2
= 4πγp−1
2 F1 (− 2 , − 2 ; 2 p; 4 ρ ) − 1 .
J ! K!

Collecting together these results yields (5.8).
Setting ρ = 1 and applying Gauss’ theorem for the value of 2 F1 (a; , b; c; 1), we deduce
e
that V(X,
X) is given by (5.9).
We remark that the absolute convergence of the series in Corollary 5.5 follows from
the absolute convergence of Gauss’ hypergeometric series. As a consequence, the series
(5.5) converges absolutely because the matrix Λ has norm less than 1.

19

Distance Correlation and Lancaster Distributions

5.3

The bivariate gamma distribution

Proposition 5.6. Suppose that the random vector (X, Y ) is distributed according to a
Sarmanov bivariate gamma distribution, as given by (3.10). Then,
V 2 (X, Y ) = 22(2−α−β)

∞
X

aj ak Aj,k (α) Aj,k (β),

(5.13)

j,k=1

where
1/2
(α)j (α)k
Aj,k (α) = 2
j! k!
Γ(2α + j + k − 1)
×
2 F1 (−j − k + 2, 1 − α − j; 2 − 2α − j − k; 2) .
Γ(α + j) Γ(α + k)
−j−k



Proof. By (3.10), there holds the expansion,
φX,Y (x, y) − φX (x) φY (y) = φX (x) φY (y)

∞
X

an L(α−1)
(x) L(β−1)
(y),
n
n

n=1

x, y > 0. Then, it follows from (4.1) that for s ∈ R,
Z ∞
Pn (s) =
exp(isx) L(α−1)
(x) φX (x) dx
n
0
Z ∞

1
exp − (1 − is)x xα−1 L(α−1)
(x) dx.
=
n
Γ(α) 0
By a direct calculation using (3.9), we obtain


(α)n
n!

1/2



(α)n
n!

1/2

Pn (s) =
=

(1 − is)−α 1 − (1 − is)−1

n

(1 − is)−(α+n) (−is)n

and, analogously,

Qn (t) =

(β)n
n!

1/2

(1 − it)−(β+n) (−it)n ,

t ∈ R.
We now calculate the integral
Z

ds
Pj (s) Pk (−s) 2 ≡
s
R



(α)j (α)k
j! k!

1/2
i

−j+k

Z
g(s)ds,

(5.14)

R

where
g(s) = sj+k−2 (1 − is)−(α+j) (1 + is)−(α+k) ,

(5.15)

20

Dueck, Edelmann, and Richards

s ∈ R. To calculate the integral on the right-hand side of (5.14), we utilize Cauchy’s
beta integral [1, p. 48] which provides that, for a, u, v ∈ C such that Re(a) > 0 and
Re(u + v) > 1,
Z
Γ(u + v − 1) u−1
(1 − is)−u (1 + ias)−v ds = 2π
a
(a + 1)1−u−v .
(5.16)
Γ(u)Γ(v)
R
To differentiate the left-hand side of (5.16) m times with respect to a, we apply the
formula,
 ∂ m
(1 + ias)−v = (−i)m sm (v)m (1 + ias)−v−m ;
∂a
by differentiating under the integral we obtain
Z
 ∂ m Z
−u
−v
m
(1 − is) (1 + ias) ds = (−i) (v)m
sm (1 − is)−u (1 + ias)−v−m ds.
∂a
R
R
To differentiate the right-hand side of (5.16) m times with respect to a, we apply
Leibniz’s formula,
m  h
 ∂ m h
i X
i
m
∂ ` u−1 i h ∂ m−`
u−1
1−u−v
a
(a + 1)
=
a
·
(a + 1)1−u−v .
∂a
`
∂a
∂a
`=0
Noting that
 
m
(−1)` (−m)`
,
=
`!
`
 ∂ `
au−1 = (−1)` au−1−` (1 − u)` ,
∂a
and
 ∂ m−`
(u + v − 1)m
,
(a + 1)1−u−v = (−1)m (a + 1)1−u−v−m+`
∂a
(2 − u − v − m)`
we obtain
 ∂ m h
∂a

a

u−1

(a + 1)

1−u−v

i
m
X

(−m)` (1 − u)`
a−` (a + 1)`
`! (2 − u − v − m)`
`=0

a + 1
= (−1)m au−1 (a + 1)1−u−v−m (u + v − 1)m 2 F1 − m, 1 − u; 2 − u − v − m;
.
a
Comparing the derivatives of the left- and right-hand sides of (5.16), we obtain
Z
sm (1 − is)−u (1 + ias)−v−m ds
m u−1

= (−1) a

(a + 1)

1−u−v−m

(u + v − 1)m

R

= 2π (−i)m au−1 (a + 1)1−u−v−m
×

Γ(u + v − 1)
Γ(u)Γ(v)


(u + v − 1)m
a + 1
F
−
m,
1
−
u;
2
−
u
−
v
−
m;
.
2 1
(v)m
a

Distance Correlation and Lancaster Distributions

21

Substituting a = 1, m = j + k − 2, u = α + j, and v = α + k − m ≡ α − j + 2, the
latter equation reduces to
Z
Γ(2α + 1)
g(s)ds = 2−2α−j−k+2 π (−i)j+k−2
Γ(α + j)Γ(α − j + 2)
R

(2α + 1)j+k−2
×
2 F1 − j − k + 2, 1 − α − j; 2 − 2α − j − k; 2 .
(α − j + 2)j+k−2
Therefore,

1/2
Z
ds
(α)j (α)k
Γ(2α + j + k − 1)
−2α+1
j−1
Pj (s) Pk (−s) 2 = 2
π (−1)
s
j! k!
Γ(α + j)Γ(α + k)
R

× 2 F1 − j − k + 2, 1 − α − j; 2 − 2α − j − k; 2 ,
and similarly for Y . Substituting these expressions into Theorem 4.1 and simplifying
the outcome, we obtain the series (5.13) as a formal expression for V 2 (X, Y ).
Finally, we verify that (5.13) converges absolutely. By (5.15),
Z
Z
|g(s)| ds =
|s|j+k−2 (1 + s2 )−(2α+j+k)/2 ds.
R

R
2

Making the change-of-variables s = t/(1 − t), the latter integral is transformed to
Z 1

1
1
t 2 (j+k−3) (1 − t)α− 2 dt = B 21 (j + k − 1), α + 21 ,
(5.17)
0

where B(·, ·) is the classical beta function, and this integral converges absolutely because
j + k − 1 > 0 and α + 1/2 > 0 for all j, k ∈ N and α > 0. Hence, to establish that
(5.13) converges absolutely, we have only to show that the series
1/2 
1/2

∞ X
∞
X
(α)k (β)k
(α)j (β)j
aj ak
(j!)2
(k!)2
j=1 k=1


× B 12 (j + k − 1), α + 21 B 21 (j + k − 1), β + 21
(5.18)
converges absolutely.
For j + k ≥ 3, it follows from (5.17) that
Z 1

1
1
1
B 2 (j + k − 1), α + 2 ≤
(1 − t)α− 2 dt =
0

1
.
α + 12

Therefore, (5.18) is bounded above by
β 2 λ2 B( 12 , α + 12 ) B( 21 , β + 12 ) +

1
(α +

1
)(β
2

+

X (β)j (β)k
λj+k
(j!) (k!)

1
)
2 j, k≥1
j+k≥3

∞

2 2

≤β λ

B( 12 , α

+

1
) B( 12 , β
2

+

1
)
2

∞

 X (β)  X (β)

1
j j
k k
+
λ
λ
(k!)
(α + 12 )(β + 12 ) j=0 (j!)
k=0

= β 2 λ2 B( 21 , α + 12 ) B( 12 , β + 21 ) +

1
(α +

1
)(β
2

+ 12 )

(1 − λ)−2β ,

22

Dueck, Edelmann, and Richards

for λ ∈ [0, 1). Hence (5.18), and then also (5.13), converges absolutely for all α, β and
all λ ∈ [0, 1).
In calculating the distance variances V(X, X) and V(Y, Y ), it is only the marginal
distributions of X and Y which are relevant. Therefore, we may assume that X and
Y have any joint distribution for which the marginal distributions are gamma with
parameters α and β, respectively.
Holding α fixed and setting β = α, the Sarmanov bivariate gamma distribution
reduces to the Kibble-Moran distribution, and the characteristic function of (X, Y ) is

−α
ψX,Y (t1 , t2 ) = (1 − it1 )(1 − it2 ) + λt1 t2
,

(5.19)

see [11, p. 436]. Next, we let λ → 1−; then ψX,Y (t1 , t2 ) converges to

−α

1 − i(t1 + t2 )
≡ E exp i(t1 + t2 )X ,
proving that if λ = 1 then X = Y , almost surely. Therefore, the distance variance
V(X, X) is a limiting case of V(X, Y ), viz.,
Z
1
ds dt
2
V (X, X) = 2
|ψX (s + t) − ψX (s)ψX (t)|2 2 2
γ1 R2
s t
Z
1
ds dt
= lim 2
|ψX,Y (s, t) − ψX (s)ψY (t)|2 2 2
λ→1− γ1 R2
s t β=α
= lim V 2 (X, Y )
λ→1−

.
β=α

Analogously, by holding β fixed and then setting α = β, we obtain
V 2 (Y, Y ) = lim V 2 (X, Y )
λ→1−

.
α=β

As a remark on the gamma distributions, note that if we replace aj in (5.13) by
[(β)j /(α)j ]1/2 λj then the result above generalizes to the distribution functions introduced by Griffiths [8]; see also (3.12) and (3.13).
As noted in Subsection 5.3, if α 6= β then Corr(X, Y ) 6≡ λ, and then it is impractical
to compare Corr(X, Y ) with ρ, the correlation coefficient in the bivariate normal case.
If α = β then Corr(X, Y ) = λ, so we will consider only the case in which α = β.
In Figure 1, we graph the difference between the distance correlation coefficient of
the Kibble-Moran bivariate gamma distribution, with shape parameter α, and the distance correlation coefficient of the bivariate normal distribution. The graphs are given
for the cases in which α = 0.1, 1, 10. Figure 1 suggests that the distance correlation
coefficient R(X, Y ) converges, as α → ∞, to the distance correlation coefficient for the
bivariate normal distribution.

23

Distance Correlation and Lancaster Distributions

Figure 1: Graphs of the difference between the distance correlation coefficient
of the Kibble-Moran bivariate gamma distribution, with shape parameter α,
and the distance correlation coefficient of the bivariate normal distribution.
The graphs are given for the cases in which α = 0.1, 1, and 10.
This result can be proved as follows: Let (X, Y )α denote a bivariate gamma random
variable with probability density function (3.10) where α = β. It follows from the characteristic function (5.19) that if (X1 , Y1 )α , . . . , (Xn , Yn )α are independent, identically
distributed random vectors with the same distribution as (X, Y )α then (X1 , Y1 )α +
· · · + (Xn , Yn )α has the same distribution as (X, Y )nα . Since (X, Y )α has finite mean
and covariance matrix then, by the Central Limit Theorem, (X, Y )nα − E(X, Y )nα converges, as n → ∞, to a bivariate normal distribution, so R((X, Y )nα ) converges to the
distance correlation coefficient of the bivariate normal distribution.
Equivalently, as α → ∞, R((X, Y )α ) converges to the distance correlation coefficient
of the bivariate normal distribution. Moreover, as the graphs indicates, the rate of
convergence is rapid. Indeed, for α = 0.1, we observe from Figure 1 that the maximum
absolute difference between R((X, Y )α ) and the distance correlation of the bivariate
normal distribution is less than 0.02; and at α = 10, the maximum absolute difference
already is negligible.

5.4

The bivariate Poisson distribution

Proposition 5.7. Suppose that the random vector (X, Y ) is distributed according to a
bivariate Poisson distribution, as given by (3.15). Then
∞
X
4j+k−1
V (X, Y ) =
(λa)j+k A2jk ,
j!
k!
j,k=1
2

(5.20)

24

Dueck, Edelmann, and Richards

where
Ajk

1
=
(j − 1)!

b(j−k)/2c 

X
`=0


j−k
(−1)` ( 21 )` ( 12 )j−`−1 1 F1 (j − ` − 21 ; j; −4a)
2`

for j ≥ k, and Ajk = Akj for j < k.
Proof. By (3.15) and (4.1), we have
Pn (s) = Qn (s) = E exp(isX) Cn (X; a),
s ∈ R. Substituting the definition (3.14) of the Poisson-Charlier polynomials Cn into
the expectation and reversing the order of summation, we obtain
Pn (s) = Qn (s) =

∞
X

exp(isx)Cn (x; a)

x=0

=

 an 1/2
n!

e −a ax
x!


(1 − e i s )n exp −a(1 − e i s ) .

Therefore, for j, k ≥ 1,
Z
ds
Pj (s)Pk (−s) 2
s
R
 aj+k 1/2 Z
 ds
=
(1 − e i s )j (1 − e −i s )k exp −a(1 − e i s + 1 − e −i s ) 2
j! k!
s
R
 aj+k 1/2 Z
 ds
(1 − e i s )j (1 − e −i s )k exp − 2a(1 − cos s) 2 .
(5.21)
=
j! k!
s
R
Changing variables in this integral from s to −s shows that the integral is symmetric
in j and k; therefore we assume, without loss of generality, that j ≥ k. We now write
(1 − e i s )j (1 − e −i s )k = (1 − e i s )j−k (1 − e i s )k (1 − e −i s )k
= (1 − e i s )j−k (2(1 − cos s))k ,
and apply the binomial theorem in the form,
(1 − e i s )j−k = (1 − cos s − i sin s)j−k

j−k 
X
j−k
=
(−i sin s)` (1 − cos s)j−k−` .
`
`=0
Then, it follows that the integral in (5.21) equals
k

2


j−k 
X
j−k
`=0

`

(−i)

`

Z

 ds
(sin s)` (1 − cos s)j−` exp − 2a(1 − cos s) 2 .
s
R

(5.22)

Distance Correlation and Lancaster Distributions

25

Expanding the exponential term,
∞
X
(−2a)m
exp − 2a(1 − cos s) =
(1 − cos s)m ,
m!
m=0



applying the half-angle identities, sin s = 2 sin 21 s cos 12 s and 1 − cos s = 2(sin 21 s)2 , and
integrating term-by-term, we deduce that (5.22) equals
k

2


j−k 
X
j−k
`=0

Z
∞
X
(−2a)m
ds
(−i)
(2 sin 21 s cos 12 s)` (2(sin 21 s)2 )j−`+m 2
`
m!
s
R
m=0

Z
j−k 
∞
X
X
(−a)m j+k+2m
j−k
ds
`
(−i)
2
(cos 21 s)` (sin 12 s)2(j+m)−` 2 . (5.23)
=
`
m!
s
R
m=0
`=0
`

If ` is odd then the latter integral is an odd function of s, so the integral equals 0.
Hence, (5.23) equals
b(j−k)/2c 


Z
∞
X
ds
j−k
(−a)m j+k+2m
2`
2
(cos 21 s)2` (sin 21 s)2(j+m−`) 2 ,
(−i)
m!
s
2`
R
m=0

X
`=0

where b(j − k)/2c denotes the greatest integer less than or equal to (j − k)/2.
Next, we introduce the formula
Z
(2` − 1)!! (2k − 3)!! π
ds
(cos 12 s)2` (sin 21 s)2k 2 =
,
(5.24)
s
(2` + 2k − 2)!! 2
R
` = 0, 1, 2 . . ., k = 1, 2, 3, . . .. This result is well-known for the case in which ` = 0 (see
[7, p. 483, 3.821(10)]), and the general case can be established by induction on ` with
the inductive step being obtained by writing
(cos 21 s)2(`+1) ≡ (cos 12 s)2` (cos 12 s)2 = (cos 12 s)2` (1 − sin2 12 s).
Hence, we find that (5.23) equals
π
2

b(j−k)/2c 

X
`=0


∞
X
j−k
(−a)m j+k+2m (2(j + m − `) − 3)!!
`
(−1) (2` − 1)!!
2
.
2`
m!
(2(j
+
m)
−
2)!!
m=0

Writing each double factorial in terms of rising factorials, and simplifying the resulting
expressions, we find that this sum equals
π
2

b(j−k)/2c 

X
`=0


∞
1
X
j−k
(−a)m 2m (j − ` − 21 )m
`
j+k−` ( 2 )j−`−1
(−1) (2` − 1)!! 2
2
2`
(j − 1)! m=0 m!
(j)m

b(j−k)/2c 
2j+k−1 X
j−k
=π
(−1)` ( 21 )` ( 21 )j−`−1 1 F1 (j − ` − 12 ; j; −4a).
(j − 1)! `=0
2`

26

Dueck, Edelmann, and Richards

Substituting this result into (5.23), we obtain
Z
ds
Pj (s)Pk (−s) 2
s
R
 aj+k 1/2 2j+k−1 b(j−k)/2c
X j − k 
(−1)` ( 21 )` ( 21 )j−`−1 1 F1 (j − ` − 12 ; j; −4a).
=π
2`
j! k!
(j − 1)! `=0
Substituting this result into Theorem 4.1 and simplifying the outcome, we obtain the
series (5.20) as a formal expression for V 2 (X, Y ).
Finally, we establish the absolute convergence of the series (5.20). On applying to
(5.21) the identity
|1 − e i s | = |1 − e −i s | = 2(1 − cos s)

1/2

= 2 sin2 21 s

1/2

and the inequality

exp − 2(1 − cos s) ≤ 1,
s ∈ R, we obtain
Z
 aj+k 1/2 Z
 ds
ds
Pj (s)Pk (−s) 2 ≤
|1 − e i s |j |1 − e i s |k exp − 2a(1 − cos s) 2
s
j! k!
s
R
R
 (4a)j+k 1/2 Z
(j+k)/2 ds
≤
.
sin2 12 s
j! k!
s2
R
By the Cauchy-Schwarz inequality,
Z
Z
(j+k)/2 ds
j/2
k/2 ds
2 1
2 1
2 1
≡
sin 2 s
sin
s
sin
s
2
2
s2
s2
R
R
Z
1/2 Z

j ds
k ds 1/2
2 1
2 1
≤
sin 2 s
sin 2 s
.
s2
s2
R
R
Since (2k − 3)!!/(2k − 2)!! ≤ 1 for all k ∈ N then it follows from (5.24) with ` = 0 that
Z
j ds
≤ π;
sin2 21 s
s2
R
therefore,
Z
Pj (s)Pk (−s)
R

 (4a)j+k 1/2
ds
≤
π,
s2
j! k!

and the same holds for the functions Qj . Substituting these bounds into the general
series expansion (4.3), we obtain the upper bound
2

V (X, Y ) ≤

∞ X
∞
X
(4λa)j+k
j=1 k=1

j! k!

2
= exp(4λa) − 1 < ∞,

27

Distance Correlation and Lancaster Distributions
for all λ ∈ [0, 1] and a > 0. Therefore, the series (5.20) converges absolutely.

To calculate the distance variance, the argument given in the bivariate gamma case
remains valid here. By Koudou [14, p. 103], the characteristic function of (X, Y ) is


ψX,Y (s, t) = exp a(1 − λ)(e i s − 1) + a(1 − λ)(e i t − 1) + aλ(e i (s+t) − 1) .
Therefore,


lim ψX,Y (s, t) = exp a(e i (s+t) − 1 ≡ ψX (s + t),

λ→1−

so we obtain
V 2 (X, X) = V 2 (Y, Y ) = lim V 2 (X, Y ).
λ→1−

5.5

The bivariate negative binomial distribution

Proposition 5.8. Suppose that the random vector (X, Y ) is distributed according to a
bivariate negative binomial distribution, as given by (3.17). Then,
∞
X
(β)j (β)k
V (X, Y ) = (1 − c)
(1 + c2 )−2β−2j 22k (λc)j+k A2jk ,
j!
k!
j,k=1
2

4β

(5.25)

where
Ajk




`
j−k 
∞
X
X
j−k
j−k
2c
(β + j) − `
`1
`2
=
(−c) (−1) (|`1 − `2 |)!
`1
`2
`!
1 + c2
` ,` =0
`=0
1

2

|`1 −`2 |

×

X
m=0

(−2)m (m)|`1 −`2 | 2k+m−1 ( 21 )k+m−1
1
2 F1 (−`, k + m − 2 ; k + m; 2),
(|`1 − `2 | − m)! (2m)! (k + m − 1)!

for j ≥ k, and Ajk = Akj for j < k.
Proof. By (3.17),
φX,Y (x, y) − φX (x) φY (y) = φX (x) φY (y)

∞
X

λn Mnβ,c (x) Mnβ,c) (y),

n=1

x, y ∈ N0 . Then, it follows from (4.1) that for s ∈ R,
Pn (s) = Qn (s) = E exp(isX) Mnβ,c (X)
∞
X
cx (β)x
=
exp(isx) Mnβ,c (x) (1 − c)β
.
x!
x=0
Substituting the definition of the Meixner polynomial as given in (3.16), we obtain
Pn (s) = (1 − c)

β

∞
 cn (β) 1/2 X
n

n!

n
k
cx (β)x X (−n)k (−x)k
exp(isx)
1 − c−1 .
x!
(β)k k!
x=0
k=0

28

Dueck, Edelmann, and Richards

Letting u ≡ ce i s and interchanging the order of summation, we obtain
∞
n
 cn (β) 1/2 X
X
ux (−x)k (β)x
(−n)k
n
−1 k
Pn (s) = (1 − c)
(1 − c )
.
n!
(β)
x!
k k!
x=0
k=0
β

(5.26)

Noting that
x

k k

∞
X
ux (−x)k (β)x

k k

k



d
du



d
du

k X
∞

d
du

k

(−x)k u = (−1) u

ux ,

we obtain

x=0

x!

= (−1) u

= (−1)k uk



x=0

ux (β)x
x!

(1 − u)−β

= (−1)k (β)k uk (1 − u)−β−k .
Substituting this result into the inner summation in (5.26) and simplifying the resulting
expression, we obtain

k
n
 n
1/2
X
(−n)k (1 − c)e i s
β c (β)n
i s −β
Pn (s) = (1 − c)
(1 − ce )
n!
k!
1 − ce i s
k=0

n
 n
1/2
(1 − c)e i s
β c (β)n
i s −β
= (1 − c)
(1 − ce )
1−
n!
1 − ce i s
 cn (β) 1/2
n
= (1 − c)β
(1 − ce i s )−β−n (1 − e i s )n .
n!
Therefore, for j, k ≥ 1,
Z
 j+k (β) (β) 1/2
ds
j
k
2β c
Pj (s)Pk (−s) 2 = (1 − c)
s
j! k!
R
Z
ds
× (1 − ce i s )−β−j (1 − ce −i s )−β−k (1 − e i s )j (1 − e −i s )k 2 .
s
R
(5.27)
Changing the variable of integration from s to −s shows that (5.27) is symmetric in j
and k; so we assume, without loss of generality, that j ≥ k.
Next, we write the integrand in (5.27) in the form
(1 − ce i s )−β−j (1 − ce −i s )−β−k (1 − e i s )j (1 − e −i s )k
= [(1 − ce i s )(1 − ce −i s )]−β−j [(1 − e i s )(1 − e −i s )]k (1 − ce −i s )j−k (1 − e i s )j−k
= (1 + c2 − 2c cos s)−β−j [2(1 − cos s)]k (1 − ce −i s )j−k (1 − e i s )j−k

−β−j
2c
k
2 −β−j
= 2 (1 + c )
1−
cos s
(1 − cos s)k (1 − ce −i s )j−k (1 − e i s )j−k .
2
1+c

Distance Correlation and Lancaster Distributions

29

By the binomial theorem,
(1 − ce

−i s j−k

)

i s j−k

(1 − e )



j−k 
X
j−k
j−k
(−c)`1 (−1)`2 e i (`2 −`1 )s ,
=
`1
`2
` ,` =0
1

2

so it follows that the integral in (5.27) is a linear combination of integrals of the form
Z 
−β−j
2c
k −i (`1 −`2 )s ds
1−
cos
s
(1
−
cos
s)
e
1 + c2
s2
R
which, after symmetrizing with s replaced by −s, equals
Z 
−β−j
1
2c
ds
1−
cos
s
(1 − cos s)k [e −i (`1 −`2 )s + e i (`1 −`2 )s ] 2
2
2 R
1+c
s
Z 

−β−j
 ds
2c
k
1−
=
cos
s
(1
−
cos
s)
cos
(`
−
`
)s
.
1
2
1 + c2
s2
R
To calculate this integral, we write
∞
−β−j X
`

(β + j)`  2c
2c
cos
s
=
cos
s
,
1−
2
1 + c2
`!
1
+
c
`=0
and apply the Chebyshev polynomials Tn , given by
|`1 −`2 |



cos (`1 −`2 )s = T`1 −`2 (cos s) = (|`1 −`2 |)!

X
m=0

(−2)m

(m)|`1 −`2 |
(1−cos s)m ;
(|`1 − `2 | − m)! (2m)!

see [7, p. 1056, 8.942(1)]. Then, we see that we need to calculate
Z
ds
(cos s)` (1 − cos s)k+m 2 .
s
R
Using the standard half-angle transformations for the cosine function and applying
(5.24), we obtain
Z
Z
ds
`
k+m ds
(cos s) (1 − cos s)
= (1 − 2 sin2 21 s)` (2 sin2 12 s)k+m 2
2
s
s
R
R


Z
`
X `
ds
= 2k+m
(−2)r (sin2 21 s)r+k+m 2
s
r
R
r=0
`  
X
`
(2r + 2k + 2m − 3)!!
k+m−1
=2
π
(−2)r
. (5.28)
r
(2r
+
2k
+
2m
−
2)!!
r=0
Expressing these double factorials in terms of rising factorials, we find that (5.28) equals
π

2k+m−1 ( 21 )k+m−1
1
2 F1 (−`, k + m − 2 ; k + m; 2).
(k + m − 1)!

Collecting together all terms, we obtain (5.25).
Finally, a proof of the absolute convergence of (5.25) can be obtained using arguments similar to those used to establish convergence in the previous subsections.

30

6

Dueck, Edelmann, and Richards

Summary and conclusions

We have derived a general series expansion for the distance covariance and distance
correlation coefficients for the class of Lancaster distributions. This result resolves the
fundamental obstacle arising in calculating the singular integrals used to define distance
correlation. We have established the utility of the result by applying it to derive
the distance correlation for the bivariate normal distribution and its generalizations
of Lancaster type, the multivariate normal distributions, and the bivariate gamma,
Poisson, and negative binomial distributions which are of Lancaster type.
In computing any of the series obtained in this paper, we can derive upper bounds on
the maximum discrepancy arising from the use of a finite number of terms of the series
by applying well-known methods of Kotz, et al. [12] together with classical bounds for
the various hypergeometric series appearing in the expansions.

Acknowledgments
The research of Dueck and Edelmann was supported by the Deutsche Forschungsgemeinschaft (German Research Foundation) within the programme “Spatio/Temporal
Graphical Models and Applications in Image Analysis,” grant GRK 1653. The research of Richards was partially supported by the U.S. National Science Foundation,
grant DMS-1309808; by a sabbatical leave-of-absence at Heidelberg University; and
by a Romberg Guest Professorship at the Heidelberg University Graduate School for
Mathematical and Computational Methods in the Sciences, funded by the German
Universities Excellence Initiative grant GSC 220/2.

References
[1] Andrews, G. E., Askey, R., and Roy, R. (1999). Special Functions. Cambridge
University Press, New York.
[2] Bar-Lev, S. K., Bshouty, D., Letac, G., Lu, I., and Richards, D. St. P. (1993).
The diagonal multivariate natural exponential families and their classification. J.
Theoret. Probab., 7, 883–929.
[3] Diaconis, P., Khare, K., and Saloff-Coste, L. (2008). Gibbs sampling, exponential
families and orthogonal polynomials. Statist. Sci., 23, 151–178.
[4] Dueck, J., Edelmann, D., Gneiting, T., and Richards, D. St. P. (2012). The affinely
invariant distance correlation. Preprint, http://arxiv.org/abs/1210.2482v1.
[5] Dueck, J., Edelmann, D., Gneiting, T., and Richards, D. St. P. (2014). The affinely
invariant distance correlation. Bernoulli, 20, 2305–2330.

Distance Correlation and Lancaster Distributions

31

[6] Dueck, J., Edelmann, D., and Richards, D. (2015). A generalization of an integral
arising in the theory of distance correlation. Statist. & Probab. Lett., 97, 116–119.
[7] Gradshteyn, I. S., and Ryzhik, I. M. (1994). Tables of Integrals, Series, and Products, fifth edition (A. Jeffrey, Editor). Academic Press, New York.
[8] Griffiths, R. C. (1969). The canonical correlation coefficients of bivariate gamma
distributions. Ann. Math. Statist., 40, 1401-1408.
[9] Huo, X., and Székely, G. J. (2016). Fast computing for distance covariance. Technometrics, 58, 435–447.
[10] Kong, J., Klein, B. E. K., Klein, R., and Wahba, G. (2012). Using distance correlation and SS-ANOVA to access associations of familial relationships, lifestyle
factors, diseases, and mortality. Proc. Natl. Acad. Sci. USA, 109, 20352–20357.
[11] Kotz, S., Balakrishnan, N., and Johnson, N. L. (2000). Continuous Multivariate
Distributions, Volume 1, Models and Applications, second edition. Wiley, New
York.
[12] Kotz, S., Johnson, N. L., and Boyd, D. W. (1967). Series representations of distributions of quadratic forms in normal variables. I. Central case. Ann.Math. Statist.,
38, 823–837.
[13] Koudou, A. E. (1996). Probabilités de Lancaster. Expo. Math., 14, 247–275.
[14] Koudou, A. E. (1998). Lancaster bivariate probability distributions with Poisson,
negative binomial and gamma margins. Test, 7, 95–110.
[15] Lancaster, H. O. (1958). The structure of bivariate distributions. Ann. Math.
Statist., 29, 719–736.
[16] Lancaster, H. O. (1969). The Chi-Squared Distribution. Wiley, New York.
[17] Lyons, R. (2013). Distance covariance in metric spaces. Ann. Probab., 41, 3284–
3305.
[18] Martı́nez-Gómez, E., Richards, M. T., and Richards, D. St. P. (2014). Distance
correlation methods for discovering associations in large astrophysical databases.
Astrophys. J., 781, 39 (11 pp.).
[19] Pommeret, D. (2004). A characterization of Lancaster probabilities with margins
in a multivariate additive class. Sankhyā, 66, 1–19.

32

Dueck, Edelmann, and Richards

[20] Richards, M. T., Richards, D. St. P., and Martı́nez-Gómez, E. (2014). Interpreting
the distance correlation results for the COMBO-17 survey. Astrophys. J. Lett.,
784, L34 (5 pp.).
[21] Rizzo, M. L., and Székely, G. J. (2010). DISCO analysis: A nonparametric extension of analysis of variance. Ann. Appl. Statist., 4, 1034–1055.
[22] Sarmanov, I. O. (1970). The approximate computation of the coefficient of correlation between functions of dependent random variables. Math. Notes Acad. Sciences
USSR, 7, 373–377.
[23] Sarmanov, O. V. (1966). Generalized normal correlation and two-dimensional
Fréchet classes. Soviet Math. Dokl., 7, 596–599.
[24] Sarmanov, O. V. (1970). A gamma-correlation process and its properties. Dokl.
Akad. Nauk SSSR, 191, 30–32.
[25] Sarmanov, O. V., and Bratoeva, Z. N. (1967). Probabilistic properties of bilinear
expansions in Hermite polynomials. Theor. Probab. Appl., 12, 470–481.
[26] Sejdinovic, D., Sriperumbudur, B., Gretton, A., and Fukumizu, K. (2013). Equivalence of distance-based and RKHS-based statistics in hypothesis testing. Ann.
Statist., 41, 2263–2291.
[27] Székely, G. J., and Rizzo, M. (2009). Brownian distance covariance. Ann. Appl.
Statist., 3, 1236–1265.
[28] Székely, G. J., and Rizzo, M. (2012). On the uniqueness of distance correlation.
Statist. & Probab. Lett., 82, 2278–2282.
[29] Székely, G. J., and Rizzo, M. (2013). The distance correlation t-test of independence in high dimension. J. Multivariate Anal., 117, 193–213.
[30] Székely, G. J., and Rizzo, M. (2014). Partial distance correlation with methods for
dissimilarities. Ann. Statist., 42, 2382–2412.
[31] Székely, G. J., Rizzo, M. L., and Bakirov, N. K. (2007). Measuring and testing
dependence by correlation of distances. Ann. Statist., 35, 2769–2794.
[32] Withers, C. S., and Nadarajah, S. (2010). Expansions for the multivariate normal.
J. Multivariate Anal., 101, 1311–1316.
[33] Zhou, Z. (2012). Measuring nonlinear dependence in time-series, a distance correlation approach. J. Time Series Anal., 33, 438–457.

