Ensemble Estimation of Mutual Information
Kevin R. Moon1 , Kumar Sricharan2 , and Alfred O. Hero III3,*

arXiv:1701.08083v1 [] 27 Jan 2017

1

Department of Genetics, Yale University — 2 Xerox PARC — 3 EECS Department, University of Michigan

Abstract—We derive the mean squared error convergence rates
of kernel density-based plug-in estimators of mutual information
measures between two multidimensional random variables X and
Y for two cases: 1) X and Y are both continuous; 2) X is continuous and Y is discrete. Using the derived rates, we propose an
ensemble estimator of these information measures for the second
case by taking a weighted sum of the plug-in estimators with
varied bandwidths. The resulting ensemble estimator achieves the
1/N parametric convergence rate when the conditional densities
of the continuous variables are sufficiently smooth. To the best of
our knowledge, this is the first nonparametric mutual information
estimator known to achieve the parametric convergence rate for
this case, which frequently arises in applications (e.g. variable
selection in classification). The estimator is simple to implement
as it uses the solution to an offline convex optimization problem
and simple plug-in estimators. A central limit theorem is also
derived for the ensemble estimator. Ensemble estimators that
achieve the parametric rate are also derived for the first case (X
and Y are both continuous) and another case 3) X and Y may
have any mixture of discrete and continuous components.

I. I NTRODUCTION
Mutual information (MI) estimation has many applications
in machine learning including MI has been used in fMRI data
processing [1], structure learning [2], independent subspace
analysis [3], forest density estimation [4], clustering [5], neuron classification [6], and intrinsically motivated reinforcement
learning [7], [8]. Another particularly common application
is feature selection or extraction where features are chosen
to maximize the MI between the chosen features X and
the outcome variables Y [9]–[12]. In many of these applications, the predictor labels have discrete components (e.g.
classification labels) while the input variables have continuous
components. To the best of our knowledge, there are currently
no nonparametric MI estimators that are known to achieve the
parametric mean squared error (MSE) convergence rate 1/N
when X and/or Y contain discrete components. Also, while
many nonparametric estimators of MI exist, most can only
be applied to specific information measures (e.g. Shannon or
Rényi information). In this paper, we provide a framework
for nonparametric estimation of a large class of MI measures
where we only have available a finite population of i.i.d.
samples. We separately consider three cases: 1) X and Y
are both continuous; 2) X is continuous and Y is discrete; 3)
X and Y may have any mixture of discrete and continuous
components. We focus primarily on the second case which
includes the problem of feature selection in classification. We
derive a MI estimator for this case that achieves the parametric
∗ This research was partially supported by ARO MURI grant W911NF-151-0479 and DOE grant de-na0002534.

MSE rate when the conditional densities of the continuous
variables are sufficiently smooth. We also show how these
estimators are extended to the first and third cases.
Our estimation method applies to other MI measures in
addition to Shannon information, which have been the focus
of much interest. The authors of [9] defined an information measure based on a quadratic divergence that could be
estimated more efficiently than Shannon information. A MI
measure based on the Pearson divergence was considered in
[13] for computational efficiency and numerical stability. The
authors of [14] and [3] used minimal spanning tree generalized
nearest-neighbor graph approaches, respectively, to estimate
Rényi information.
A. Related Work
Many estimators for Shannon MI between continuous random variables have been developed. A popular k-nn-based
estimator was proposed in [15] which is a modification of
the entropy estimator derived in [16]. However, these estimators only achieve the parametric convergence rate when
the dimension of each of the random variables is less than 3
[17]. Similarly, the Rényi information estimator in [3] does not
achieve the parametric rate. Some other estimators are based
on maximum likelihood estimation of the likelihood ratio [18]
and minimal spanning trees [19].
Recent work has focused on nonparametric divergence estimation for purely continuous random variables. One approach
[20]–[23] uses an optimal kernel density estimator (KDE) to
achieve the parametric convergence rate when the densities are
at least d [22], [23] or d/2 [20], [21] times differentiable where
d is the dimension of the data. These optimal KDEs require
knowledge of the density support boundary and are difficult
to construct near the boundary. Numerical integration may
also be required for estimating some divergence functionals
under this approach, which can be computationally expensive.
In contrast, our approach to MI estimation does not require numerical integration and can be performed without knowledge
of the support boundary.
More closely related work [24]–[28] uses an ensemble
approach to estimate entropy or divergence functionals. These
works construct an ensemble of simple plug-in estimators by
varying the neighborhood size of the density estimators. They
then take a weighted average of the estimators where the
weights are chosen to decrease the bias with only a small
increase in the variance. The parametric rate of convergence is
achieved when the densities are either d [24]–[26] or (d+1)/2
[27], [28] times differentiable. These approaches are simple to

implement as they only require simple plug-in estimates and
the solution of an offline convex optimization problem. These
estimators have also performed well in various applications
[29]–[33]
Finally, the authors of [34] showed that k-nn or KDE based
approaches underestimate the MI when the MI is large. As
MI increases, the dependencies between random variables increase which results in less smooth densities. Thus a common
approach to overcome this issue is to require the densities to
be smooth [20]–[28].
B. Contributions
In the context of this related work, we make the following
novel contributions in this paper: (1) For continuous random
variables (case 1), we extend the asymptotic bias and variance
results for divergence estimators [27], [28] to kernel density
plug-in MI estimators without boundary correction [35] by
incorporating machinery to handle the dependence between
the product of marginal density estimators (Section II), (2) we
extend the theory to handle discrete random variables in the
mixed cases (cases 2 and 3) by reformulating the densities
as a mixture of the conditional density of the continuous
variables given the discrete variables (Section III), and (3) we
leverage this theory for the mixed cases in conjunction with
the generalized theory of ensemble estimators [27], [28] to
derive, to the best of our knowledge, the first non-parametric
estimator that achieves a parametric rate of MSE convergence
of O (1/N ) for the mixed cases (Section IV), where N is
the number of samples available from each distribution. We
also derive a central limit theorem for the ensemble estimators
(Section IV-C). We verify the theory through experiments
(Section V).
II. C ONTINUOUS R ANDOM VARIABLES
In this section, we obtain MSE convergence rates of plugin MI estimators when X and Y are continuous (case 1 in
Section I). This will enable us to derive the MSE convergence
rates of plug-in MI estimators when X is continuous and
Y is discrete and when X and Y may have any mixture
of continuous and discrete components (respectively, cases
2 and 3 in Section I). These rates can then be used to
derive ensemble estimators that achieve the parametric MSE
rate. Let fX (x), fY (y), and fXY (x, y) be dX , dY , and

dX + dY = d-dimensional densities. Let g(t1 , t2 ) = g tt12
(e.g. g(t1 , t2 ) = log(t1 /t2 ) for Shannon information). We
define a family of MIs as

Z 
fX (x)fY (y)
fXY (x, y)dxdy.
(1)
G1 (X; Y) = g
fXY (x, y)

and KY (·) be kernel functions with ||KX ||∞ , ||KY ||∞ < ∞
where ||K||∞ = supx |K(x)|. The KDE for fX is


N
Xj − Xi
1 X
.
(2)
KX
f̃X,hX (Xj ) =
hX
M hdXX i=1
i6=j

The KDEs f̃Y,hY (Yj ) and f̃Z,hZ (Xj , Yj ) (where hZ =
(hX , hY )) for estimating fY and fXY , respectively, are defined similarly using KY and the product kernel KX · KY .
Then G1 (X; Y) is estimated as
!
N
1 X
f̃X,hX (Xi )f̃Y,hY (Yi )
G̃hX ,hY =
.
(3)
g
N i=1
f̃Z,hZ (Xi , Yi )
B. Convergence Rates

To derive the convergence rates of G̃hX ,hY we assume that
1) fX , fY , fXY , and g are smooth; 2) fX and fY have
bounded support sets SX and SY ; 3) fX , fY , and fXY are
strictly lower bounded on their support sets. More specifically,
we assume that the densities belong to the bounded Hölder
class Σ(s, H) (the precise definition is included in the appendices) which implies that the densities are r = ⌊s⌋ times
differentiable. These assumptions are comparable to those
in similar studies on asymptotic convergence analysis [20]–
[26], [28]. To derive the convergence rates without boundary
corrections, we also assume that 4) the boundary of the support
set is smooth with respect to the corresponding kernels. The
full assumptions are
• (A.0): The kernels KX and KY are symmetric product
kernels with bounded support.
• (A.1): There exist constants ǫ0 , ǫ∞ such that 0 < ǫ0 ≤
fX (x) ≤ ǫ∞ < ∞, ∀x ∈ SX , ǫ0 ≤ fY (y) ≤ ǫ∞ , ∀y ∈
SY , and ǫ0 ≤ fXY (x, y) ≤ ǫ∞ , ∀(x, y) ∈ SX × SY .
• (A.2): Each of the densities belong to Σ(s, H) in the
interior of their support sets with s ≥ 2.
• (A.3): g (t1 /t2 ) has an infinite number of mixed derivatives wrt t1 and t2 .
∂ k+l g(t1 ,t2 )
/(k!l!), k, l = 0, 1, . . . are strictly
• (A.4):
l
∂tk
1 ∂t2
upper bounded for ǫ0 ≤ t1 , t2 ≤ ǫ∞ .
• (A.5): Let K be either KX or KY , S either SX or SY , h
either hX or hY . Let px (u) : Rd → R be a polynomial in
u of order q ≤ r = ⌊s⌋ whose coefficients are a function
of x and are r − q times differentiable. For any positive
integer t
!t
Z
Z
K(u)px (u)du

x∈S

where vt (h) admits the expansion

A. The KDE Plug-in Estimator
When both X and Y are continuous with marginal densities
fX and fY , the MI functional G1 (X; Y) can be estimated
using KDEs. Assume that N i.i.d. samples {Z1 , . . . , ZN } are
available from the joint density fXY with Zi = (Xi , Yi )T . Let
M = N − 1 and let hX , hY be kernel bandwidths. Let KX (·)

dx = vt (h),

u:K(u)>0, x+uh∈S
/

vt (h) =

r−q
X
i=1


ei,q,t hi + o hr−q ,

for some constants ei,q,t .
Assumption (A.5) states that the support of the density is
smooth with respect to the kernel K in the sense that the

expectation with respect to any random variable u of the area
of the kernel that falls outside the support S is a smooth
function of the bandwidth h provided that the distribution
function px (u) of u is smooth (e.g. s ≥ 2). The inner integral
captures this expectation while the outer integral averages
this inner integral over all points near the boundary of the
support. The vt (h) term captures the fact that the smoothness
of this expectation is proportional to the smoothness of the
function px (u). As an example, this smoothness assumption
is satisfied when the support is rectangular and the kernel is the
uniform rectangular kernel [27], [28]. Note that this boundary
assumption does not result in parametric convergence rates for
the plug-in estimator G̃hX ,hY , which is in contrast with the
boundary assumptions in [20]–[23]. However, the estimators in
[20]–[23] perform boundary correction, which requires knowledge of the density support boundary and complex calculations
at the boundary in addition to the boundary assumptions, to
achieve the parametric convergence rates. In contrast, we use
ensemble methods to improve the resulting convergence rates
of G̃hX ,hY without boundary correction.
Theorem 1. (Bias) Under assumptions A.0 − A.5 and for
general g, the bias of G̃hX ,hY is
r
r
X
X

h
i
B G̃hX ,hY
=

c10,i,j hiX hjY +

j=0 i=0
i+j6=0

+O hsX + hsY +

c11
N hdXX hdYY

1
N hdXX hdYY

If g (t1 , t2 ) also has j, l-th order mixed derivatives
β
tα
1 t2

!

.

∂ j+l
∂tj1 ∂tl2

(4)
that

depend on t1 and t2 only through
for some α, β ∈ R for
each 1 ≤ k, l, ≤ λ, the bias of G̃hX ,hY is
h
i
B G̃hX ,hY
⌊λ/2⌋

=

X

r
X

m,n=0
i,j=0
i+j+m+n6=0

+

⌊λ/2⌋ r
r
X XX
m=1 i=0 j=0

c11,j,i,m,n 

hiX hjY
m 
n
N hdXX
N hdYY


m
c13,m,n,j hiX hjY / N hdXX hdYY



λ/2 
+O hsX + hsY + 1/ N hdXX hdYY
.

(5)

The constants in both (4) and (5) depend only on the densities and their derivatives, the functional g and its derivatives,
and the kernels. They are independent of N, hX , and hY .
The purpose of Theorem 1 is two-fold. First, we use
Theorem 1 to derive the bias expressions for the MI plugin estimators when X and Y may have a mixture of discrete
and continuous components (cases 2 and 3) in Section III.
Second, in conjunction with Theorem 2 which follows, the
results in Theorem 1 can be used to derive MI ensemble
estimators in Appendix B-A that achieve the parametric MSE
convergence rate when the densities are sufficiently smooth.

The expression in (5) enables us to achieve the parametric rate
under less restrictive smoothness assumptions on the densities
(s > d/2 for (5) compared to s ≥ d for (4)). The extra
condition required on the mixed derivatives of g to obtain the
expression in (5) is satisfied, for example, for Shannon and
Renyi information measures.
Theorem 2. (Variance) If the functional g is Lipschitz continuous in both of its arguments with Lipschitz constant Cg ,
then the variance of G̃hX ,hY is
h
i 22C 2 ||KX · KY ||2
g
∞
.
V G̃hX ,hY ≤
N
Similar to Theorem 1, Theorem 2 is used to derive variance
expressions for the MI plug-in estimators under cases 2 and
3. Theorem 2 is also necessary to derive optimally weighted
ensemble estimators. The proofs of Theorems 1 and 2 are
similar to the proofs of the bias and variance results for
the divergence functional estimators in [27]. The primary
difference is in handling certain products of the marginal
KDEs that appear in the expansion of the MSE. See Appendix
C and D for details.
Theorems 1 and 2 indicate that for the MSE of the plug-in
estimator to go to zero for case 1, we require hX , hY → 0
and N hdXX hdYY → ∞. The Lipschitz assumption on g is
comparable to other nonparametric estimators of distributional
functionals [20]–[23], [27]. Specifically, assumption A.1 ensures that functionals such as those for Shannon and Renyi
informations are Lipschitz on the space ǫ0 to ǫ∞ .
III. M IXED R ANDOM VARIABLES
In this section, we extend the results of Section II to MI
estimation when X and Y may have a mixture of discrete and
continuous components. For simplicity, we focus primarily on
the important case when X is continuous and Y is discrete
(case 2 in Section I). The more general case when X and Y
may have any mixture of continuous and discrete components
(case 3 in Section I) is discussed in Section III-B. As an
example of the former case, if Y is a predictor variable (e.g.
classification labels), then the MI between X and Y indicates
the value of X as a predictor of Y. Although Y is discrete,
fXY = fZ is also a density. Let SX be the support of the
density fX and SY be the support of the probability mass
function fY . The MI is
G2 (X; Y)
X Z  fX (x)fY (y) 
fXY (x, y)dx
(6)
=
g
fXY (x, y)
y∈SY

Z 
X
fX (x)
fX|Y (x|y)dx.
=
fY (y) g
fX|Y (x|y)
y∈SY
PN
Let Ny = i=1 1{Yi =y} where y ∈ SY . Let f̃X,hX be as
in (2) and define Xy = {Xi ∈ {X1 , . . . , XN } |Yi = y}. Then
if Xi ∈ Xy , the KDE of fX|Y (x|y) is


X
1
Xi − Xj
f̃X|y,hX|y (Xi ) =
.
KX
X
hX|y
(Ny − 1) hdX|y
Xj ∈Xy
i6=j

We define the plug-in estimator G̃hX ,hX|Y of (6) as

1 X 
g f̃X,hX (X)/f̃X|y,hX|y (X) ,
G̃hX ,hX|y =
Ny
X∈Xy

X Ny
=
G̃hX ,hX|y .
N

=⇒ G̃hX ,hX|Y

(7)

y∈SY

A. Convergence Rates
To apply the theory of optimally weighted ensemble estimation to G̃hX ,hX|Y , we need to know its MSE as a function
of the bandwidths and the sample size.
Theorem 3. (Bias) Assume that assumptions A.0 − A.5 apply
to the functional g, the kernel KX , and the densities fX and
fX|Y . Assume that hX|y = lN−β
with 0 < β < d1X and l a
y
positive number. Then the bias of G̃hX ,hX|Y is
i
h
B G̃hX ,hX|Y
=

r
r
X
X

c13,i,j hiX lj N −jβ +

j=0 i=0
i+j6=0

+ O hsX + N −sβ

c14,X
c14,Y
+ dX 1−βdX
l N
N hdXX

1
1
+ 1−βdX
+
dX
N
N hX

!

.

If g (t1 , t2 ) also has j, l-th order mixed derivatives
β
tα
1 t2

depend on t1 and t2 only through
each 1 ≤ j, l ≤ λ, then the bias is
h
i
B G̃hX ,hX|Y
⌊λ/2⌋

=

that

for some α, β ∈ R for

hiX lj N −jβ
m
c14,j,i,m,n 
n
N hdXX
(ldX N 1−βdX )
m,n=0
i,j=0
i+j+m+n6=0



1
1


+ O hsX + N −sβ + 
.
λ/2 +
λ/2
1−βd
X)
(N
N hdXX

(9)

Proof: We focus on (8) as (9) follows similarly. It can be
shown that


i
i
h
X Ny h
B G̃hX ,hX|y Y1 , . . . , YN  .
B G̃hX ,hX|Y = E 
N
y∈SY

The conditional bias of G̃hX ,hX|y given Y1 , . . . , YN can then
be obtained from Theorem 1 as
i
h
B G̃hX ,hX|y Y1 , . . . , YN
=

r
X

c10,i,j hiX hjX|y

i,j=0
i+j6=0

+O

hsX

+

hsX|y

1
1
+
+
dX
X
N y hX
Ny hdX|y

!

i=0

(10)

From [36], the i-th central moment of Ny has the form of
i ⌊i/2⌋
h
X
E (NY − N fY (y))i =
cn,i (fY (y))N n .
n=0

has terms proportional to N 1−γ−i+n ≤
N
for i = 0, 1, . . . since n ≤ ⌊i/2⌋. Then since
there is an N in the denominator of (7), this leaves terms of
the form of N −γ when i = 0, 1 and N −1 for i ≥ 2. This
completes the proof for the bias. See Appendix E for more
details.
Thus E N1−γ
y
1−γ−⌊i/2⌋

r
X

X

α

Nα
y = (Ny − N fY (y) + N fY (y))

∞ 
X
α
α−i
i
=
(N fY (y))
(Ny − N fY (y)) ,
i
i=0
 
=⇒ E Nα
y

∞ 
i
h
X
α
α−i
i
(N fY (y))
E (Ny − N fY (y)) .
=
i



(8)
∂ j+l
∂tj1 ∂tl2

Then given that hX|y ∝ N−β
y , (7) gives terms of the form
of N1−γ
with
γ
>
0.
N
is a binomial random variable
y
y
with parameter fY (y), N trials, and mean N fY (y). Thus we
need to compute the fractional moments of a binomial random
variable. By the generalized binomial theorem, we have that



Theorem 4. If the functional g is Lipschitz continuous in
both of its arguments and SY is finite, then the variance of
G̃hX ,hX|Y is O(1/N ).
Proof: By the law of total variance, we have
ii
i
h h
h
V G̃hX ,hX|Y = E V G̃hX ,hX|Y Y1 , . . . , YN
ii
h h
+ V E G̃hX ,hX|Y Y1 , . . . , YN .

Given all of the Yi ’s, the estimators G̃hX ,hX|y are all independent since they use differenthsets of Xi ’s for each y. From
i
Theorem 2, we know that V G̃hX ,hX|Y Y1 , . . . , YN =
P

2
N
/N
O
. Taking the expectation then yields
y
y∈SY
O(1/N ).
For the
h second term, we know ifrom the proof of Theorem 3
that E G̃hX ,hX|Y Y1 , . . . , YN yields a sum of terms of
the form of Nγy /N for 0 < γ ≤ 1. Taking the variance
of the sum
 of these terms yields a sum of terms of the
form V Nγy /N 2 (the covariance terms can be bounded by
the Cauchy-Schwarz
inequality to yield similar terms). Then

V Nγy can be bounded by taking a Taylor series expansion of
the functions Nγy and N2γ
y at the point N fY (y) which yields
an expression that depends
 the central moments of Ny .
 on
From this, we obtain V Nγy = O(N ) which completes the
proof. See Appendix E for details.
Theorems 3 and 4 provide exact expressions for the bias and
bounds on the variance of the plug-in MI estimator, respectively. It is shown in Section IV that the MSE of the plug-in
estimator converges very slowly to zero under this setting.

However, Theorems 3 and 4 provide with us the necessary
information for applying the theory of optimally weighted
ensemble estimation to obtain estimators with improved rates.
This is done in Section IV.
B. Extension to Other Cases
The results in Section III-A can be extended to the case
where X and/or Y may have a mixture of continuous and
discrete components (case 3 in Section I). This scenario can be
divided further into three different cases: A) X is continuous
and Y has a mixture of discrete and continuous components;
B) X and Y both have a mixture of discrete and continuous
components; C) Y is discrete and X has a mixture of discrete
and continuous components. Consider case A first. Denote the
discrete and continuous components of Y as Y1 and Y2 ,
respectively. Denote the respective support sets as SY1 and
SY2 . We can then write
G3A (X; Y)
X Z  fX (x)fY (y1 , y2 ) 
fXY (x, y1 , y2 )dxdy2
=
g
fXY (x, y1 , y2 )
y1 ∈SY1

Z 
X
fX (x)fY2 |Y1 (y2 |y1 )
=
fY1 (y1 ) g
fXY2 |Y1 (x, y2 |y1 )
y1 ∈SY1

× fXY2 |Y1 (x, y2 |y1 )dxdy2 .

(11)

The subscript 3A indicates that we are considering case
A under the third case described in the introduction. The
expression in (11) is very similar to the expression in (6).
After plugging in KDEs for the corresponding densities and
conditional densities, a nearly identical procedure to that in
Section III-A can be followed to derive the bias and variance
of the corresponding plug-in estimator.
Now consider case B. Denote the discrete and continuous
components of X as X1 and X2 , respectively. Then if Y1 is
the discrete component of Y, then the expression inside the
g functional in (11) includes fX1 (x1 )fY1 (y1 )/fX1 Y1 (x1 , y1 ).
Thus the plug-in estimator must include estimators for
fX1 (x1 ), fY1 (y1 ), and fX1 Y1 (x1 , y1 ). Define Ny1 =
P
N
i=1 1{Y1,i =y1 } where Y1,i is the discrete component of Yi .
Then the estimator we use for fY1 (y1 ) is Ny1 /N . The estimators for fX1 (x1 ) and fX1 Y1 (x1 , y1 ) are defined similarly.
The bias and variance expressions of this plug-in estimator can
then be derived with some slight modifications of Theorems 1
and 2. See Appendix E-C for an expression for G3B (X; Y) in
this case and a sketch of these modifications. Case C follows
similarly as the expression inside the g functional in (11)
includes fX1 (x1 )fY (y)/fX1 Y (x1 , y) where all the terms are
probability mass functions.
The resulting bias and variance expressions in these settings
are analogous to those in Theorems 1, 2, and 3 as the variance
will be O(1/N ) and the bias will depend on expansions of the
bandwidths for the various KDEs. Ensemble methods can then
be applied to improve the MSE convergence rates as described
in the next section.

IV. E NSEMBLE E STIMATION OF MI
A. Mixed Random Variables
We again focus on the case where X is continuous and
Y is discrete (case 2 in Section I). If no bias correction
is performed, then Theorem 3 shows that the optimal bias

rate of the plug-in estimator G̃hX ,hX|Y is O 1/N 1/(dX +1) ,
which converges very slowly to zero when dX is not small.
We use the theory of optimally weighted ensemble estimation
to improve this rate. An ensemble of estimators is formed by
choosing different bandwidth values. Consider first the case
where (8) applies. Let L be a set of real positive numbers with
|L| = L. This set will parameterize the bandwidths for f̃X,hX
and f̃X|y,hX|y resulting in L estimators in the ensemble. While
different parameter sets for f̃X,hX and f̃X|y,hX|y can be chosen,
we only use one set here for simplicity
√ of exposition. To ensure
that the final terms in (8) are O(1/ N ) when s ≥ d, for each
estimator in the ensemble we choose hX (l) = lN −1/(2dX )
−1/(2dX )
and hX|y (l) = lNy
where l ∈ L. Define
P w to be a
weight vector parameterized by l ∈ L with l∈L w(l) = 1
and define
X
X Ny
G̃hX (l),hX|y (l) .
(12)
G̃w,1 =
w(l)
N
l∈L

y∈SY

From Theorem 3, the bias of G̃w,1 is
r


i XX
h
−i
θ w(l)li N 2dX
B G̃w,1 =
l∈L i=1

+O

√
 −s

−1
L||w||2 N 2dX + N 2
,

(13)

where we use θ notation to omit the constants.
We use the general theory of optimally weighted ensemble
estimation in [28] to improve the MSE convergence rate of
the plug-in estimator by using the weights to cancel
o lower
n the
order terms in (13). The theory is as follows. Let Êl
be
l∈L
an indexed ensemble
of
estimators
with
the
weighted
ensemble
P
estimator Êw = l∈L w(l)Êl satisfying:
•

C.1. Let ci be constants depending on the underlying
density, J = {i1 , . . . iI } a finite index set with I < L,
ψi (l) basis functions depending only on the parameter l
and not on N , φi (N ) functions of the sample size N that
are independent of l. Assume the bias is


h i X
1
.
ci ψi (l)φi (N ) + O √
B Êl =
N
i∈J

•

C.2. Assume the variance is
 
 
h i
1
1
+o
.
V Êl = cv
N
N

Theorem 5. [28] If conditions
C.1 and C.2 hold for an
n o
ensemble of estimators Êl
, then there exists a weight
l∈L

vector w0 such that the MSE of Êw0 attains the parametric

rate of convergence of O (1/N ). The weight w0 is the solution
to the offline convex optimization problem
minw ||w||
P 2
subject to
l∈L w(l)
P = 1,
γw (i) = l∈L w(l)ψi (l) = 0, i ∈ J.

(14)

To apply Theorem 5 to an ensemble of estimators,
√ all φi (N )
functions that converge to zero slower than 1/ N and the
corresponding ψi (l) functions must be known for the base
estimator. Otherwise, Theorem 5 can only be guaranteed to
improve the bias up to the slowest unknown bias rate. This
theorem was applied in [28] to the problem of divergence
functional estimation where the plug-in estimator has slowly
converging bias but the resulting ensemble estimator achieves
the parametric rate for sufficiently smooth densities.
We apply Theorem 5 to the ensemble estimator G̃w,1 as
conditions C.1 and C.2 are satisfied with φi (N ) = N −i/(2dX )
and ψi (l) = li for i ∈ {1, . . . r} as seen in (8) and (13). If
s ≥ dX , then the MSE of the optimally weighted estimator
G̃w0 ,1 is O(1/N ). A similar approach can be used for the
case where X contains a mixture of continuous and discrete
components and Y is discrete (or vice versa). To the best of
our knowledge, these are the first nonparametric estimators
to achieve the MSE parametric rate in this setting of mixed
random variables.
If the mixed derivatives of the functional g satisfy the
extra condition required for (9), we can define an ensemble
estimator G̃w0 ,2 that achieves the parametric MSE rate if
s > dX /2. For simplicity, we focus primarily on G̃w0 ,1 . See
Appendix B-B for details on G̃w0 ,2 .
In practice, the optimization problem in (14) typically
results in a very large increase in variance. Thus we follow
the lead of [24]–[27] and use a relaxed version of (14):
minw ǫP
subject to
l∈L w(l) = 1,
1

γw (i)N 2 φi (N ) ≤ ǫ, i ∈ J,

(15)

2

kwk2 ≤ η.

As shown in [24]–[27], the ensemble estimator G̃w0 ,1 using
the resulting weight vector from the optimization problem in
(15) still achieves the parametric MSE convergence rate under
the same assumptions as described previously. It was also
shown in [27] that the heuristic of setting η = ǫ works well
in practice. Algorithm 1 summarizes the estimator G̃w0 ,1 .
A similar approach can be used to derive an ensemble
estimator G̃cont
w0 ,1 for the case when X and Y are continuous
(case 1 in Section I). See Appendix B-A for details. The
case where X and Y both contain a mixture of discrete and
continuous components follows similarly.

Algorithm 1 Optimally weighted ensemble MI estimator
G̃w0 ,1
Input: L positive real numbers L, samples {Z1 , . . . , ZN }
from fXY , dimension dX , function g, kernel KX
Output: The optimally weighted MI estimator G̃w0 ,1
1: Solve for w0 using (15) with basis functions ψi (l) = l i ,
φi (N ) = N −i/(2dX ) , l ∈ L, and 0 ≤ i ≤ dX .
2: for all l ∈ L and y ∈ SY do
P
3:
Ny ← N
i=1 1{Yi =y}
−1/(2dY )
4:
hX (l) ← lN −1/(2dX ) , hX|y (l) ← lNy
5:
for Xi ∈ Xy do
6:
Calculate f̃X,hX (l) (Xi ), f̃X|y,hX|y (l) (Xi ) as described
in the text
7:
end for


P
f̃
(X)
8:
G̃hX (l),hX|y (l) ← N1y X∈Xy g f̃ X,hX (l) (X)
X|y,hX|y (l)

9:

10:

end for P
P
G̃w0 ,1 ← l∈L w0 (l) y∈Sy

Ny
N G̃hX (l),hX|y (l)

1) Select the minimum and maximum bandwidth parameter
to produce density estimates that satisfy the following:
first the minimum bandwidth should not lead to a zerovalued density estimate at any sample point; second
the maximum bandwidth should be smaller than the
diameter of the support.
2) Ensure the bandwidths are sufficiently distinct. Similar
bandwidth values lead to negligible decrease in bias and
many bandwidth values may increase ||w0 ||2 resulting in
an increase in variance [24].
3) Select L = |L| > |J| = I to obtain a feasible solution
for the optimization problems in (14) and (15). We find
that choosing a value of 30 ≤ L ≤ 60, and setting L to
be L linearly spaced values between the minimum and
maximum values described above works well in practice.
The resulting ensemble estimators are robust in the sense that
they are not sensitive to the exact choice of the bandwidths
or the number of estimators as long as the the rough rules-ofthumb given above are followed. Moon et al [27], [28] gives
more details on ensemble estimator parameter selection for
continuous divergence estimation. These details also apply to
the continuous parts of the mixed cases for MI estimation in
this paper.
Since the optimal weight w0 can be calculated offline, the
computational complexity of the estimators is dominated by
the construction
of the KDEs which has a complexity of

O N 2 using the standard implementation. For very large
datasets, more efficient KDE implementations (e.g. [37]) can
be used to reduce the computational burden.

B. Parameter Selection
In theory, the theoretical results of the previous sections
hold for any choice of the bandwidth vectors as determined
byL. In practice, we find that the following rules-of-thumb
for tuning the parameters lead to high-quality estimates in the
finite sample regime.

C. Central Limit Theorem
We finish this section with central limit theorems for the
ensemble estimators. This enables us to perform hypothesis
testing on the mutual information.

Theorem 6. Let G̃cont
be a weighted ensemble estimator
w
when X and Y are continuous with bandwidths hX (lX ) and
hY (lY ) for each estimator in the ensemble. Assume that the
functional g is Lipschitz in both arguments with Lipschitz
constant Cg and that hX (lX ), hY (lY ) = o(1), N → ∞,
and N hdXX (lX ), N hdYY (lY ) → ∞ for each lX ∈ LX and
lY ∈ LY . Then for fixed LX and LY , and if S is a standard
normal random variable,
!
i

h
i r h

kernel and the optimally weighted ensemble estimator G̃w0 ,1
for various sample sizes and for d = 4, 6, 9, respectively. The
ensemble estimator outperforms the standard plug-in estimator,
especially for larger sample sizes and larger dimensions.
This demonstrates that while an individual kernel estimator
performs poorly, an ensemble of estimators including the
individual estimator performs well.
For the second case, Y has six possible outcomes (i.e.
|SY | = 6) and respective probabilities Pr(Y = 0) = 0.35,
≤ t → Pr (S ≤ t) . Pr(Y = 1) = 0.2, Pr(Y = 2) = Pr(Y = 3) = 0.15,
Pr
G̃cont
− E G̃cont
/ V G̃cont
w
w
w
Pr(Y = 4) = 0.1, and Pr(Y = 5) = 0.05. We chose
α = 0.5 and d = 6. The conditional covariances matrices
The proof is based on an application of Slutsky’s Theorem
are again 0.1 × Id and the conditional means are, respectively,
preceded by an application of the Efron-Stein inequality (see
µ̄0 = 0.25 × 1̄d , µ̄1 =
 0.75 × 1̄d , and µ̄2 =T 0.5 × 1̄d , µ̄T3T=
Appendix F).
T
T T
0.25
×
1̄
,
0.5
×
1̄
, µ̄4 = 0.75 × 1̄2 , 0.375 × 1̄4 ,
4
2
If the space SY is finite, then the ensemble estimators for

T
T T
the mixed component case also obey a central limit theorem. and µ̄5 = 0.5 × 1̄4 , 0.25 × 1̄2 . The parameters for the
The proof follows by an application of Slutsky’s Theorem ensemble estimators and the KDE plug-in estimators are the
same as in the top three plots in Figure 1. The bottom plot in
combined with Theorem 6.
Figure 1 again compares the ensemble estimator to the plug-in
Corollary 7. Let G̃w be a weighted ensemble estimator when KDE estimator. The ensemble estimator also outperforms the
X is continuous and Y is discrete with bandwidths hX (l) plug-in estimator in this setting.
and hX|y (l) for each estimator in the ensemble. Assume that
the functional g is Lipschitz in both arguments and that
VI. C ONCLUSION
X
hX , hX|y = o(1), N → ∞, and N hdXX , N hdX|y
→ ∞ for
We derived the MSE convergence rates for plug-in KDEeach l ∈ L and ∀y ∈ SY with Sy finite. Then for fixed L,
based estimators of MI measures between X and Y when
!
they have only continuous components and for the case where

i r h
h
i
Pr

G̃w − E G̃w

/ V G̃w ≤ t

→ Pr (S ≤ t) .

V. E XPERIMENTAL VALIDATION
In this section, we validate our theory by estimating the
Rényi-α MI integral (i.e. g(x) = xα in (6); see [38]) where X
is a mixture of truncated Gaussian random variables restricted
to the unit cube and Y is a categorical random variable. We
choose Rényi MI as it has received recent interest (e.g. [3]) and
the estimation problem does not reduce to entropy estimation
in contrast with Shannon MI. Thus this is a clear case where
there are no other nonparametric estimators that are known to
achieve the parametric MSE rate.
We consider two cases. In the first case, Y has three
possible outcomes (i.e. |SY | = 3) and respective probabilities
Pr(Y = 0) = Pr(Y = 1) = 2/5 and Pr(Y = 2) = 1/5.
The conditional covariance matrices are all 0.1 × Id and
the conditional means are, respectively, µ̄0 = 0.25 × 1̄d ,
µ̄1 = 0.75 × 1̄d , and µ̄2 = 0.5 × 1̄d , where Id is the d × d
identity matrix and 1̄d is a d-dimensional vector of ones.
This experiment can be viewed as the problem of estimating
MI (e.g. for feature selection or Bayes error bounds) of a
classification problem where each discrete value corresponds
to a distinct class, the distribution of each class overlaps
slightly with others, and the class probabilities are unequal. We
use α = 0.5. We set L to be 40 linearly spaced values between
1.2 and 3. The bandwidth in the KDE plug-in estimator is also
set to 2.1N −1/(2d).
The top three plots in Figure 1 shows the MSE (200 trials) of
the plug-in KDE estimator of the MI integral using a uniform

Y is discrete and X is continuous. We also showed how
convergence rates can be obtained for the case when X and/or
Y contain a mixture of discrete and continuous components.
Using these rates, we defined ensemble estimators that achieve
an MSE rate of O(1/N ) when the densities are sufficiently
smooth and showed that a central limit theorem also holds. To
the best of our knowledge, this is the first nonparametric MI
estimator that achieves the MSE convergence rate of O(1/N )
in this setting of mixed random variables (i.e. X and Y are
not both purely discrete or purely continuous).
R EFERENCES
[1] B. Chai, D. Walther, D. Beck, and L. Fei-Fei, “Exploring functional
connectivities of the human brain using multivariate information analysis,” in Advances in neural information processing systems, pp. 270–278,
2009.
[2] K. R. Moon, M. Noshad, S. Y. Sekeh, and A. O. Hero III, “Information
theoretic structure learning with confidence,” in Proc IEEE Int Conf
Acoust Speech Signal Process, 2017.
[3] D. Pál, B. Póczos, and C. Szepesvári, “Estimation of rényi entropy and
mutual information based on generalized nearest-neighbor graphs,” in
Advances in Neural Information Processing Systems, pp. 1849–1857,
2010.
[4] H. Liu, L. Wasserman, and J. D. Lafferty, “Exponential concentration for
mutual information estimation with application to forests,” in Advances
in Neural Information Processing Systems, pp. 2537–2545, 2012.
[5] J. Lewi, R. Butera, and L. Paninski, “Real-time adaptive informationtheoretic optimization of neurophysiology experiments,” in Advances in
Neural Information Processing Systems, pp. 857–864, 2006.
[6] E. Schneidman, W. Bialek, and M. J. B. II, “An information theoretic
approach to the functional classification of neurons,” Advances in Neural
Information Processing Systems, vol. 15, pp. 197–204, 2003.
[7] S. Mohamed and D. J. Rezende, “Variational information maximisation
for intrinsically motivated reinforcement learning,” in Advances in
Neural Information Processing Systems, pp. 2116–2124, 2015.

|S |=3, d=4
Y

Kernel
Weighted

MSE

10 -2

10 -3
10 2

10 3

Sample Size N
|S |=3, d=6
Y

MSE

10 -1

Kernel
Weighted

10 -2

10 -3
10 2

10 3

Sample Size N
|SY|=3, d=9

MSE

10 -1

Kernel
Weighted

10 -2

10 -3
10 2

10 3

Sample Size N
|SY|=6, d=6

MSE

10 -1

10 -2

10 -3
10 2

Kernel
Weighted

10 3

Sample Size N
Figure 1. MSE log-log plots as a function of sample size for the uniform
kernel plug-in MI estimator ("Kernel") and the proposed optimally weighted
ensemble estimator G̃w0 ,1 ("Weighted") for the distributions described in the
text. The top three plots each correspond to the first case where |SY | = 3
and the bottom plot corresponds to the second case where |SY | = 6. The
ensemble estimator outperforms the kernel plug-in estimator, especially for
larger sample sizes. Note also that as the dimension increases, the performance
gap between the two estimators increases.

[8] C. Salge, C. Glackin, and D. Polani, “Changing the environment based
on empowerment as intrinsic motivation,” Entropy, vol. 16, no. 5,
pp. 2789–2819, 2014.
[9] K. Torkkola, “Feature extraction by non parametric mutual information
maximization,” The Journal of Machine Learning Research, vol. 3,
pp. 1415–1438, 2003.
[10] J. R. Vergara and P. A. Estévez, “A review of feature selection methods
based on mutual information,” Neural Computing and Applications,
vol. 24, no. 1, pp. 175–186, 2014.
[11] H. Peng, F. Long, and C. Ding, “Feature selection based on mutual information criteria of max-dependency, max-relevance, and minredundancy,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 27, no. 8, pp. 1226–1238, 2005.
[12] N. Kwak and C.-H. Choi, “Input feature selection by mutual information
based on parzen window,” Pattern Analysis and Machine Intelligence,
IEEE Transactions on, vol. 24, no. 12, pp. 1667–1671, 2002.
[13] M. Sugiyama, “Machine learning with squared-loss mutual information,”
Entropy, vol. 15, no. 1, pp. 80–112, 2012.
[14] J. A. Costa and A. O. Hero, “Geodesic entropic graphs for dimension
and entropy estimation in manifold learning,” IEEE Transactions on
Signal Processing, vol. 52, no. 8, pp. 2210–2221, 2004.
[15] A. Kraskov, H. Stögbauer, and P. Grassberger, “Estimating mutual
information,” Physical review E, vol. 69, no. 6, p. 066138, 2004.
[16] L. Kozachenko and N. N. Leonenko, “Sample estimate of the entropy
of a random vector,” Problemy Peredachi Informatsii, vol. 23, no. 2,
pp. 9–16, 1987.
[17] W. Gao, S. Oh, and P. Viswanath, “Demystifying fixed k-nearest neighbor information estimators,” arXiv preprint arXiv:1604.03006, 2016.
[18] T. Suzuki, M. Sugiyama, J. Sese, and T. Kanamori, “Approximating
mutual information by maximum likelihood density ratio estimation.,”
FSDM, vol. 4, pp. 5–20, 2008.
[19] S. Khan, S. Bandyopadhyay, A. R. Ganguly, S. Saigal, D. J. Erickson III,
V. Protopopescu, and G. Ostrouchov, “Relative performance of mutual
information estimation methods for quantifying the dependence among
short and noisy data,” Physical Review E, vol. 76, no. 2, p. 026209,
2007.
[20] A. Krishnamurthy, K. Kandasamy, B. Poczos, and L. Wasserman, “Nonparametric estimation of renyi divergence and friends,” in Proceedings of
The 31st International Conference on Machine Learning, pp. 919–927,
2014.
[21] K. Kandasamy, A. Krishnamurthy, B. Poczos, L. Wasserman, and
J. Robins, “Nonparametric von mises estimators for entropies, divergences and mutual informations,” in Advances in Neural Information
Processing Systems, pp. 397–405, 2015.
[22] S. Singh and B. Póczos, “Exponential concentration of a density
functional estimator,” in Advances in Neural Information Processing
Systems, pp. 3032–3040, 2014.
[23] S. Singh and B. Póczos, “Generalized exponential concentration inequality for rényi divergence estimation,” in Proceedings of the 31st
International Conference on Machine Learning (ICML-14), pp. 333–
341, 2014.
[24] K. Sricharan, D. Wei, and A. O. Hero, “Ensemble estimators for
multivariate entropy estimation,” Information Theory, IEEE Transactions
on, vol. 59, no. 7, pp. 4374–4388, 2013.
[25] K. R. Moon and A. O. Hero, “Ensemble estimation of multivariate
f-divergence,” in Information Theory (ISIT), 2014 IEEE International
Symposium on, pp. 356–360, IEEE, 2014.
[26] K. R. Moon and A. O. Hero, “Multivariate f-divergence estimation with
confidence,” in Advances in Neural Information Processing Systems,
pp. 2420–2428, 2014.
[27] K. R. Moon, K. Sricharan, K. Greenewald, and A. O. Hero, “Nonparametric ensemble estimation of distributional functionals,” arXiv preprint
arXiv:1601.06884v2, 2016.
[28] K. R. Moon, K. Sricharan, K. Greenewald, and A. O. Hero, “Improving
convergence of divergence functional ensemble estimators,” in 2016
IEEE International Symposium on Information Theory (ISIT), 2016.
[29] Z. Szabó and A. Lőrincz, “Distributed high dimensional information
theoretical image registration via random projections,” Digital Signal
Processing, vol. 22, no. 6, pp. 894–902, 2012.
[30] S. V. Gliske, K. R. Moon, W. C. Stacey, and A. O. Hero, “The intrinsic
value of HFO features as a biomarker of epileptic activity,” in IEEE
International Conference on Acoustics, Speech, and Signal Processing,
pp. 6290–6294, 2016.

[31] K. Moon, V. Delouille, and A. O. Hero, “Meta learning of bounds on
the Bayes classifier error,” in IEEE Signal Processing and SP Education
Workshop, pp. 13–18, IEEE, 2015.
[32] K. R. Moon, J. J. Li, V. Delouille, R. De Visscher, F. Watson, and
A. O. Hero, “Image patch analysis of sunspots and active regions. I.
Intrinsic dimension and correlation analysis,” Journal of Space Weather
and Space Climate, vol. 6, no. A2, 2016.
[33] K. R. Moon, V. Delouille, J. J. Li, R. De Visscher, F. Watson, and
A. O. Hero, “Image patch analysis of sunspots and active regions. II.
Clustering via matrix factorization,” Journal of Space Weather and Space
Climate, vol. 6, no. A3, 2016.
[34] S. Gao, G. Ver Steeg, and A. Galstyan, “Efficient estimation of mutual information for strongly dependent variables,” in Proceedings of
the Eighteenth International Conference on Artificial Intelligence and
Statistics, pp. 277–286, 2015.
[35] R. J. Karunamuni and T. Alberts, “On boundary correction in kernel
density estimation,” Statistical Methodology, vol. 2, no. 3, pp. 191–212,
2005.
[36] J. Riordan, “Moment recurrence relations for binomial, poisson and
hypergeometric frequency distributions,” The Annals of Mathematical
Statistics, vol. 8, no. 2, pp. 103–111, 1937.
[37] V. C. Raykar, R. Duraiswami, and L. H. Zhao, “Fast computation of
kernel estimators,” Journal of Computational and Graphical Statistics,
vol. 19, no. 1, pp. 205–220, 2010.
[38] J. C. Principe, Information theoretic learning: Renyi’s entropy and kernel
perspectives. Springer Science & Business Media, 2010.
[39] B. Efron and C. Stein, “The jackknife estimate of variance,” The Annals
of Statistics, pp. 586–596, 1981.
[40] R. Durrett, Probability: Theory and Examples. Cambridge University
Press, 2010.
[41] A. Gut, Probability: A Graduate Course. Springer Science & Business
Media, 2012.

A PPENDIX A
H ÖLDER C LASS
We derive MSE convergence rates for the plug-in estimators in terms of the smoothness of the densities which we characterize
by the Hölder Class.
Pd
Definition 1 (Hölder Class). Let X ⊂ Rd be a compact space. For r = (r1 , . . . , rd ), ri ∈ N, define |r| = i=1 ri and
|r|
Dr = ∂xr1∂...∂xrd . The Hölder class Σ(s, H) of functions on L2 (X ) consists of the functions f that satisfy
1

d

s−r

|Dr f (x) − Dr f (y)| ≤ H kx − yk

,

for all x, y ∈ X and for all r s.t. |r| ≤ ⌊s⌋.
For notation, let EZ denote the conditional expectation given Z.
A PPENDIX B
MI E NSEMBLE E STIMATION E XTENSIONS
A. Continuous Random Varables
We can also apply Theorem 5 to obtain MI estimators that achieve the parametric rate for the case when X and Y are
dX dY
−1/2
continuous.
for the O(1/(N hdXX hdYY )) terms to
√ For general g, (4) in the main paper indicates that we need hX hY ∝ N
be O(1/ N ). We consider the more general case where the parameters may differ for hX and hY . Let LX and LY be sets of
real, positive numbers with |LX | = LX and |LY | = LY . For each estimator in the ensemble, choose
P lX ∈ LX and lY ∈ LY and
set hX (lX ) = lX N −1/(2(dX +dY )) and hY (lY ) = lY N −1/(2(dX +dY )) . Define the matrix w s.t. lX ∈LX ,lY ∈LY w(lX , lY ) = 1.
i j
From Theorems 1 and 2, conditions C.1 and C.2 are satisfied if s ≥ dX + dY with ψi,j (lX , lY ) = lX
lY and φi,j (N ) =
−(i+j)/(2(dX +dY ))
N
for 0 ≤ i, j ≤ dX + dY s.t. 0 < i + j ≤ dX + dY . The optimal weight w0 is calculated using (14) in the
main paper. The resulting estimator
X
G̃cont
w0 (lX , lY )G̃hX (lX ),hY (lY )
w0 ,1 =
lX ∈LX ,lY ∈LY

achieves the parametric MSE rate when s ≥ dX + dY .
Again, if the mixed derivatives of the functional g also satisfy the extra condition required for (5) in the main paper, then
we can define an estimator that achieves the parametric MSE rate under less strict smoothness assumptions. See Appendix
B-B2.
B. The ODin2 Estimators
The estimators G̃w0 ,1 and G̃cont
w0 ,1 are analogous to the ODin1 estimators in [27], [28]. In this section, we derive ensemble
estimators of MI that achieve the parametric rate under less strict smoothness assumptions on the densities. These estimators
are analogous to the ODin2 estimators in [27], [28].
1) Mixed Random Variables: We first consider the case where X is continuous and Y is discrete. Recall that if hX|y = lN−β
y
j+l
with 0 < β < d1X and l a positive number, and if g (t1 , t2 ) has j, l-th order mixed derivatives ∂t∂j ∂tl that depend on t1 and t2
1

2

β
only through tα
1 t2 for some α, β ∈ R for each 1 ≤ j, l ≤ λ, then the bias of the plug-in estimator for this case is
⌊λ/2⌋

i
h
B G̃hX ,hX|Y =

r
X

hiX lj N −jβ
m
c14,j,i,m,n 
n
N hdXX
(ldX N 1−βdX )
m,n=0
i,j=0
i+j+m+n6=0


X

1
1


.
+ O hsX + N −sβ + 
λ/2 +
λ/2 
1−βd
X
dX
)
(N
N hX

(16)

−1/(dX +δ)
Choose L to be a set of real positive numbers and let δ > 0. For each estimator in the ensemble,
√ set hX (l) = lN
−1/(dX +δ)
and hX|y (l) = lNy
where l ∈ L. This ensures that the final terms in (29) are O(1/ N ) if s ≥ (dX + δ)/2 and
λ ≥ dX /δ + 1. Define G̃w,2 as in (12) in the main paper with the chosen values of hX (l) and hX|y (l). Theorem 5 can be
−i−mδ

applied in this case as conditions C.1 and C.2 are satisfied with φi,m (N ) = N dX +δ and ψi,m (l) = li−mdX for i ∈ {0, . . . , r},
m ∈ {0, . . . ⌊λ/2⌋}, and 0 ≤ i + mδ ≤ (dX + δ)/2. Then if s ≥ (dX + δ)/2 and λ ≥ dX /δ + 1, the MSE of the optimally
weighted estimator G̃w0 ,2 is O(1/N ). Then since δ can be chosen arbitrarily close to zero, the parametric rate can be achieved
theoretically as long as s > dX /2.

The analogous divergence functional estimators for G̃w0 ,1 and G̃w0 ,2 in [27], [28] were referred to as the ODin1 and
ODin2 estimators, respectively, where ODin stands for Optimally Weighted Distributional Functional estimators. The ODin2
estimator has better statistical properties as the parametric rate is guaranteed under less restrictive smoothness assumptions on
the densities. On the other hand, the number of parameters required for the optimization problem in (14) in the main paper
is larger for the ODin2 estimator than the ODin1 estimator. In theory, this could lead to larger variance although this wasn’t
necessarily true in practice according to the experiments in [27].
2) Continuous Random Variables: We now consider the case where both X and Y are continuous. Again, if g (t1 , t2 ) has
j+l
β
j, l-th order mixed derivatives ∂t∂j ∂tl that depend on t1 and t2 only through tα
1 t2 for some α, β ∈ R for each 1 ≤ k, l, ≤ λ,
1

2

then the bias of G̃hX ,hY is

h

B G̃hX ,hY

i

⌊λ/2⌋

r
X

X

=

m,n=0
i,j=0
i+j+m+n6=0
⌊λ/2⌋

+

r X
r
X X

m=1 i=0 j=0

c11,j,i,m,n 

hiX hjY
m 
n
N hdXX
N hdYY

m

c13,m,n,j hiX hjY / N hdXX hdYY



λ/2 
dX dY
s
s
+O hX + hY + 1/ N hX hY
.

(17)

Set δ > 0 and choose hX (lX ) = lX N −1/(dX +dY +δ) and hY (lY ) = lY N −1/(dX +dY +δ) . Then conditions C.1 and C.2 are
i−mdX j−ndY
satisfied if s ≥ (dX + dY + δ)/2 and λ ≥ (dX + dY + δ)/δ with ψ1,i,j,m,n (lX , lY ) = lX
and φ1,i,j,m,n (N ) =
lY
N

−

i+j+m(dY +δ)+n(dX +δ)
dX +dY +δ

dX +dY +δ
2
dX +dY +δ
.
The
2

for 0 < i + j + m(dY + δ) + n(dX + δ) ≤

− d i+j+mδ
X +dY +δ

i−mdX j−mdY
and the terms ψ2,i,j,m (lX , lY ) = lX
lY

optimal weight w0 is again calculated using (14)
and φ2,i,j,m (N ) = N
for m ≥ 1 and i + j + mδ ≤
in the main paper and the resulting estimator G̃cont
achieves
the
parametric
MSE
convergence rate when s ≥ (dX + dY + δ)/2.
w0 ,2
Since δ can be chosen arbitrarily close to zero, the parametric rate can be achieved theoretically as long as s > (dX + dY )/2.
G̃cont
w0 ,2 is the ODin2 estimator for continuous random variables.
P ROOF

A PPENDIX C
T HEOREM 1 (B IAS )

OF

The proof of the bias results in Theorem 1 share some similarities with the proof of the bias results for the divergence
functional estimators in in [27]. The primary differences deal with the product of the marginal KDEs that appear in the
expansion of the bias terms.
The bias of G̃hX ,hY can be expressed as
!
"

#
h
i
fX (X)fY (Y)
f̃X,hX (X)f̃Y,hY (Y)
−g
B G̃hX ,hY
= E g
fXY (X, Y)
f̃Z,hZ (X, Y)
h
h
i 
i


!
(Y)
f̃
(X)
E
f̃
E
Y,h
Y
X,h
X
Y
X
f̃X,hX (X)f̃Y,hY (Y)

= E g
−g
f̃Z,hZ (X, Y)
EX,Y f̃Z,hZ (X, Y)
h
h
i
i

 


EX f̃X,hX (X) EY f̃Y,hY (Y)
f
(X)f
(Y)
X
Y
,
−g
(18)
+E g 
fXY (X, Y)
EX,Y f̃Z,hZ (X, Y)
where X and Y are drawn jointly from fXY . We can view these terms as a variance-like component (the first term) and a
bias-like component, where the respective Taylor series expansions depend on variance-like or bias-like
terms of the KDEs.

EX [f̃X,hX (X)]EY [f̃Y,hY (Y)]
We first consider the bias-like term, i.e. the second term in (18). The Taylor series expansion of g
f̃
(X,Y)
E
X,Y Z,hZ

around fX (X)fY (Y) and fXY (X, Y) gives an expansion with terms of the form of
h
h
i
i
h
i
i

EX f̃X,hX (X) EY f̃Y,hY (Y) − fX (X)fY (Y) ,
BiZ f̃X,hX (X)f̃Y,hY (Y) =
h
i
i

BiZ f̃Z,hZ (X, Y) =
EX,Y f̃Z,hZ (X, Y) − fXY (X, Y) .

(19)

Since we are not doing boundary correction, we need to consider separately the cases when Z is in the interior of the support
SX × SY and when Z is close to the boundary of the support. For precise definitions, a point Z = (X, Y ) ∈ SX × SY is



′ 
′ 
′
Y −Y
K
= 0, and a point Z ∈ SX × SY is near the
in the interior of SX × SY if for all Z ∈
/ SX × SY , KX X−X
Y
hX
hY
boundary of the support if it is not in the interior.
It can be shown by Taylor series expansions of the probability densities that for Z = (X, Y) drawn from fXY in the interior
of SX × SY , then
⌊s/2⌋

h
i
EX f̃X,hX (X)

=

h
i
EY f̃Y,hY (Y)

=

i
h
EX,Y f̃Z,hZ (Z)

=

fX (X) +

X

s
cX,j (X)h2j
X + O (hX ) ,

X

s
cY,j (Y)h2j
Y + O (hY ) ,

(20)

j=1

⌊s/2⌋

fY (Y) +

j=1

⌊s/2⌋ ⌊s/2⌋

fXY (X, Y) +

X X

2j
s
s
cXY,i,j (X, Y)h2i
X hY + O (hX + hY ) .

i=0 j=0
i+j6=0

For a point near the boundary of the support, we extend the expectation beyond the support of the density. As an example
if X is near the boundary of SX , then we get


Z
h
i
1
X−V
EX f̃i,hi (X) − fi (X) =
fX (V )dV − fX (X)
K
X
hX
hdXX V :V ∈SX
"
#


Z
1
X−V
=
KX
fX (V )dV − fX (X)


hX
hdXX V :KX X−V
>0
hX
#
"


Z
X−V
1
fX (V )dV
− dX
KX
hX
hX V :V ∈S
/ X
=

T1,X (X) − T2,X (X).

(21)

We only evaulate the density fX and its derivatives at points within the support when we take its Taylor series expansion.
Thus the exact manner in which we define the extension of fX does not matter as long as the Taylor series remains the same
and as long as the extension is smooth. Thus the expected value of T1,X (X) gives an expression of the form of (20). For the
T2,X (X) term, we can use multi-index notation on the expansion of fX to show that
#
"


Z
X−V
1
fX (V )dV
KX
T2,X (X) =
hX
hdXX V :V ∈S
/ X
Z
KX (u)fX (X + hX u)du
=
u:hX u+X∈S
/ X ,KX (u)>0

X h|α| Z
X
=
KX (u)Dα fX (X)uα du + o(hrX ).
α! u:hX u+X∈S
/ X ,KX (u)>0
|α|≤r

Then since the |α|th derivative of fX is r − |α| times differentiable, we apply the condition in assumption A.5 to obtain
E [T2,X (X)] =

r
X

ei hiX + o (hrX ) .

i=1

Similar expressions can be found for f̃Y,hY and f̃Z,hZ and for when (21) is raised to a power t. Applying this result gives for
the second term in (18),
r
r
X
X
c10,i,j hiX hjY + O (hsX + hsY ) .
(22)
j=0 i=0
i+j6=0

For the first term in (18), a Taylor series expansion of g



f̃X,hX (X)f̃Y,hY (Y)
f̃Z,hZ (X,Y)



h
h
i
i
around EX f̃X,hX (X) EY f̃Y,hY (Y) and

EX,Y f̃Z,hZ (X, Y) gives an expansion with terms of the form of
h
iq

ẽqZ,hZ (Z) =
f̃Z,hZ (Z) − EZ f̃Z,hZ (Z)
,
h
h
iq
i

ẽqXY,hX ,hY (Z) =
.
f̃X,hX (X)f̃Y,hY (Y) − EX f̃X,hX (X) EY f̃Y,hY (Y)

(23)

We can take the expected value of these expressions to obtain terms of the form of
1
1
1
1
,
,
,
N hdXX N hdYY N 2 hdXX hdYY N hdXX hdYY

(24)

and their respective powers. This can be seen for ẽqXY,hX ,hY (Z) as follows. Define










Xi − X
Xi − X
Yj − Y
Yj − Y
Vi,j (Z) = KX
KY
− EX KX
EY KY
hX
h
hX
hY
h Y′
i
= ηij (Z) − EX [ηi (X)] EY ηj (Y) .
We can then write

ẽXY,hX ,hY (Z) =

N
N X
X

1
N 2 hdXX hdYY

The binomial theorem then gives

k  
h ′
ik−l
X
 l

k
.
=
(Z) EX [ηi (X)] EY ηj (Y)
EZ ηij
l

k
(Z)
EZ Vi,j



Vi,j (Z).

i=1 j=1



(25)

l=0

By using a similar Taylor series analysis as before, for Z in the interior,

⌊s/2⌋


X
 l

2dX dY
dX 2dY
2n
(Z) = hdXX hdYY
EZ ηij
cXY,2,m,n,l (Z)h2m
h
+
O
h
h
+
h
h
.
X
Y
X
Y
X Y
m,n=0

Combining this with (20) and (25) gives

⌊s/2⌋


X
 k

dX dY
2dX dY
dX 2dY
2n
EZ Vi,j (Z) = hX hY
cXY,3,m,n,k (X)h2m
,
X hY + O hX hY + hX hY

(26)

m,n=0

where the constants depend on the densities, their derivatives, and the moments of the kernels. As an example, let q = 2. Then
due to the independence between the Zi samples,


EZ ẽ2XY,hX ,hY (Z) =
=
=

N
X

1
X 2dY
N 4 h2d
X hY

1

E



Z
X 2dY
N 2 h2d
X hY
⌊s/2⌋

1

N 2 hdXX hdYY

EZ [Vi,j (Z)Vm,n (Z)]

i,j,m,n=1

X

m,n=0


2
Vi,j
(Z) +

(N − 1)

X 2dY
N 2 h2d
X hY

EZ [Vi,j (Z)Vi,n (Z)]
⌊s/2⌋

2n
cXY,3,m,n,2 (X)h2m
X hY

+

X

1
X

m,n=0 i,j=0
i+j6=0

cXY,4,m,n,i,j (X)

2n
h2m
X hY
X jdY
N hid
X hY

+O



1
N



,

where the last step follows from (26) and a similar analysis of EZ [Vi,j (Z)Vi,n (Z)]. For q > 2, it can be shown that if n(q)
is the set of integer divisors of q including 1 but excluding q, then




 

i ⌊s/2⌋
h
X
X
cXY,6,i,j,q,m,n (Z)
1
 X cXY,5,i,j,q,n (Z)
 2i 2j
q
EZ ẽXY,hX ,hY (Z) =
.
h
h
+
O
+




q−n 
q−n
q−m
X Y


N
d
d
d
d
X
X
Y
Y
i,j=0 n∈n(q) N 2 hX hY
N hY
m∈n(q)∪{q} N hX

n∈n(q)∪{q}
m+n6=2q

i
h
A similar procedure can be used to find the expression for EZ ẽqZ,hZ (Z) . When Z is near the boundary of the supposrt, we
n
can obtain similar expressions by following a similar procedure as in the derivation of (22). This results in powers of hm
X hY
2m 2n
instead of hX hY .
h
h
i
i
For general functionals g, we can only guarantee that the mixed derivatives of g evaluated at EX f̃X,hX (X) EY f̃Y,hY (Y)

and EX,Y f̃Z,hZ (X, Y) converge to the mixed derivative evaluated at fX (X)fY (Y) and fXY (X, Y) at some rate o(1). Thus
we are left with the following terms in the bias:
!
1
1
o
+
N hdXX
N hdYY

j+l

β
However, if we know that g (t1 , t2 ) has j, l-th order mixed derivatives ∂t∂j ∂tl that depend on t1 and t2 only through tα
1 t2 for
1
2
some α, β ∈ R, then by the generalized binomial theorem, we find that

m

⌊s/2⌋
∞  
α

X
X
α α−m
s 
ci,j (X)h2j
fX (X) 
.
EX f̃X,hX (X) =
X + O (hX )
m
m=0
j=1
α

α

A similar result holds for EY f̃Y,hY (Y)
and EZ f̃Z,hZ (Z) . Combining these expressions with 24 completes the proof.
P ROOF

OF

A PPENDIX D
T HEOREM 2 (VARIANCE )

As for the bias, the proof of the variance result in Theorem 2 is similar to the proof of the variance result in [27] and so we
do not present all of the details. The primary differences again deal with the product of the marginal KDEs. The proof uses
the Efron-Stein inequality [39]:
′

′

Lemma 8. (Efron-Stein Inequality) Let X1 , . . . , Xn , X1 , . . . , Xn be independent random variables on the space S. Then if
f : S × · · · × S → R, we have that
V [f (X1 , . . . , Xn )] ≤


n
2 
′
1X
E f (X1 , . . . , Xn ) − f (X1 , . . . , Xi , . . . , Xn )
.
2 i=1

In this case we consider the samples {Z1 , . . . , ZN } and

′

G̃hX ,hY . By the triangle inequality,
′

G̃hX ,hY − G̃hX ,hY

′

Z1 , Z2 . . . , ZN

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

1
g
N

≤

n

N2
1 X
g
+
N j=2

!

o

and the respective estimators G̃hX ,hY and

′

′

f̃X,hX (X1 )f̃Y,hY (Y1 )
′
′
f̃Z,hZ (X1 , Y1 )

−g

f̃X,hX (Xj )f̃Y,hY (Yj )
f̃Z,hZ (Xj , Yj )

!

′

−g

!

′

f̃X,hX (Xj )f̃Y,hY (Y1 )
′

f̃Z,hZ (X1 , Y1 )

!

.

By the Lipschitz condition on g, the first term in (27) can be decomposed into terms of the form of
′

f̃Z,hZ (Z1 ) − f̃Z,hZ (Z1 ) ,
′

′

f̃X,hX (X1 )f̃Y,hY (Y1 ) − f̃X,hX (X1 )f̃Y,hY (Y1 ) .
By making a substitution in the expectation, it can be shown that
E



′

f̃Z,hZ (Z1 ) − f̃Z,hZ (Z1 )

2



≤ 2||KX · KY ||2∞ .



X1 − Xi
hX

For the product of the marginal KDEs, we have that

f̃X,hX (X1 )f̃Y,hY (Y1 ) =
=

1
M 2 hdXX hdYY

N X
N
X
i=2 j=2

KX

1
1
f̃Z,hZ (Z1 ) +
d
2
M
M hXX hdYY




Y1 − Yj
KY
hY




X
Y1 − Yj
X1 − Xi
KY
.
KX
hX
hY
i6=j



(27)

By applying the triangle inequality, Jensen’s inequality, and similar substitutions, we get
E



′

′

f̃X,hX (X1 )f̃Y,hY (Y1 ) − f̃X,hX (X1 )f̃Y,hY (Y1 )

2




2
′
2
f̃Z,hZ (Z1 ) − f̃Z,hZ (Z1 )
E
M2
2(M − 1)
+
×
X 2dY
M 3 h2d
hY
X




X 
X1 − Xi
Y1 − Yj
E KX
KY
hX
hY
i6=j
!
!!2 
′
′
Y1 − Yj
X1 − Xi

KY
−KX
hX
hY


≤

4 + 2(M − 1)2
||KX · KY ||2 .
M2

≤
For the second term in (27), it can be shown that
E



′

f̃Z,hZ (Zi ) − f̃Z,hZ (Zi )

2








Y1 − Yj
X1 − Xi
KY
E KX
X 2dY
hX
hY
M 2 h2d
X hY
!
!!2 
′
′
Y1 − Yj
X1 − Xi

KY
−KX
hX
hY
1

=

2||KX · KY ||2∞
.
M2

≤
By a similar approach,

′

′

f̃X,hX (Xi )f̃Y,hY (Yi ) − f̃X,hX (Xi )f̃Y,hY (Yi )

=

1

′

f̃Z,hZ (Zi ) − f̃Z,hZ (Zi ) +

M 2 hdXX hdYY





+

X

KX

n=2
n6=i

=⇒ E



Xi − Xn
hX



KY



X
KY

n=2
n6=i

Yi − Y1
hY





Yi − Yn
hY

′

Yi − Y1
hY

− KY

′

′



f̃X,hX (Xi )f̃Y,hY (Yi ) − f̃X,hX (Xi )f̃Y,hY (Yi )

2



KX

!!



Xi − X1
hX



′

− KX

Xi − X1
hX




,

≤ 6||KX · KY ||2∞



1
(M − 2)2
+
2
M
M4



We can then apply the Cauchy Schwarz inequality to bound the square of the second term in (27) to get

N2
 X
E 
g
j=2

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

!


! 2
f̃X,hX (X1 )f̃Y,hY (Y1 )
2
2
 
 ≤ 14Cg ||KX · KY ||∞ .
′
f̃Z,hZ (X1 , Y1 )
′

−g

′

Applying Jensen’s inequality in conjunction with these results gives
E



′

G̃hX ,hY − G̃hX ,hY

Applying the Efron-Stein inequality finishes the proof.

2



≤

44Cg2 ||KX · KY ||2∞
.
N2

!!

T HEORY

FOR

A PPENDIX E
M IXED R ANDOM VARIABLES

A. Proof of Theorem 3 (Bias)
1
dX .

Let hX|y = lN−β
for some positive l and 0 < β <
y
bias of the plug-in estimator G̃hX ,hX|Y

i
h
=
B G̃hX ,hX|Y

r
r
X
X

Under assumptions A.0 − A.5, we prove that for general g, the

c13,i,j hiX lj N −jβ +

j=0 i=0
i+j6=0

+O

hsX

+N

Furthermore, if g (t1 , t2 ) has j, l-th order mixed derivatives
α, β ∈ R, then for any positive integer λ ≥ 2, the bias is
i
h
=
B G̃hX ,hX|Y

r
r
X
X

c13,i,j hiX lj N −jβ +

j=0 i=0
i+j6=0
λ/2

+

r
X

r
X



−sβ

+

1
N hdXX

∂ j+l
∂tj1 ∂tl2

c14,X
c14,y
+ dX 1−βdX
dX
l N
N hX

+

1
N 1−βdX

1
+
N

!

.

β
that depend on t1 and t2 only through tα
1 t2 for some

λ/2 λ/2 r
r
XXX
X

n −nβ
hm
Xl N
c14,j,i,m,n 
j
i
j=1 i=1 m=0 n=0
N hdXX (ldX N 1−βdX )


n −nβ
hm
hm ln N −nβ


Xl N
c14,m,n,j,X X
j + c14,m,n,j,Y d
j
X N 1−βdX )
dX
(l
j=1 m=0 n=0
N hX



X


+O hsX + N −sβ + 

1

N hdXX

(28)

λ/2 +

1

λ/2
(N 1−βdX )

+

1
.
N

(29)

We only prove (28) as the proof of (29) is identical. The bias of G̃hX ,hX|Y is
i
h
B G̃hX ,hX|Y

i
h
= E G̃hX ,hX|Y − G(X; Y)




X Ny
fX (X)

= E
G̃hX ,hX|y − g
N
fX|Y (X|Y)
y∈SY
 



X Ny
f
(X)
X
= E E 
Y, Y1 , . . . , YN 
G̃hX ,hX|y − g
N
fX|Y (X|Y)
y∈SY





X Ny  
f
(X)
X
Y, Y1 , . . . , YN 
E
= E
G̃hX ,hX|y − g
N
fX|Y (X|Y)
y∈SY


i
X Ny h
B G̃hX ,hX|y Y1 , . . . , YN  ,
= E
N
y∈SY

where we use the law of total expectation and the fact that

P

y∈SY

Ny
N

= 1. Let hX|y = lN−β
for some positive l and
y

0<β<

1
dX .

From Theorem 1, the conditional bias of G̃hX ,hX|y given Y1 , . . . , YN is
i
h
=
B G̃hX ,hX|y Y1 , . . . , YN

r
r
X
X

c10,i,j hiX hjX|y +

j=0 i=0
i+j6=0

c11,X
c11,y
+
dX
X
N y hX
Ny hdX|y

1
1
+
+O hsX + hsX|y +
dX
X
N y hX
Ny hdX|y
=

r
r
X
X

c10,i,j hiX lj N−jβ
+
y

j=0 i=0
i+j6=0

+O

hsX

+

N−sβ
y

!

c11,y
c11,X
+
dX
X
d
X
N y hX
l N1−βd
y

1
1
+
+ 1−βdX
Ny hdXX
Ny

!

.

(30)

Ny is a binomial random variable Multiplying (30) by Ny results in terms of the form of N1−γ
with γ ≥ 0. Ny is a binomial
y
random variable with parameter fY (y),N trials, and mean N fY (y). We can compute the fractional moments of a binomial
random variable by using the generalized binomial theorem to obtain (see the main paper)

∞ 
i
h
X
 α
α
(N fY (y))α−i E (NY − N fY (y))i
=
E Ny
i
i=0

=


⌊i/2⌋
∞ 
X
X
α
α−i
(N fY (y))
cn,i (fY (y))N n
i
n=0

i=0

=


⌊i/2⌋
∞ 
X
X
α
α−i
fY (y)
cn,i (fY (y))N α−i+n ,
i
n=0

i=0

where we use the following expression for the i-th central moment of a binomial random variable derived by Riordan [36]:
h
i ⌊i/2⌋
X
i
E (NY − N fY (y)) =
cn,i (fY (y))N n .
n=0

If α = 1 − γ, then dividing by N results in terms of the form of N −γ−i+n . Since n ≤ ⌊i/2⌋, −γ − i + n is always less than
zero and is only greater than −1 if i = 0. This completes the proof.
B. Proof of Theorem 4 (Variance)
As for the bias, we assume that hX|y = lNy−β for some positive l and 0 < β < d1X . By the law of total variance, we have
ii
ii
h h
i
h h
h
(31)
V G̃hX ,hX|Y = E V G̃hX ,hX|Y Y1 , . . . , YN + V E G̃hX ,hX|Y Y1 , . . . , YN .

Note that given all of the Yi ’s, the estimators G̃hX ,hX|y are all independent since they use different sets of Xi ’s for each y.
By Theorem 2, we have


i
h
X N2y 1

= O
V G̃hX ,hX|Y Y1 , . . . , YN
·
N 2 Ny
y∈SY


X Ny
.
= O
N2
y∈SY


1

Taking the expectation wrt Y1 , . . . YN then gives O N for the first term in (31).
For the second term in (31), from (30) we have that for general g


r
i
h
X
1
X
+ N−sβ
+ N1−βd
= O
N−jβ
+
E G̃hX ,hX|y Y1 , . . . , YN
y
y
y
N
y
j=0
= O (f (Ny )) .

′

By the Efron-Stein inequality, we have that if Ny is an independent and identically distributed realization of Ny , then



 ′ 2 
X Ny
′
1 X


f (Ny )
≤
V
E Ny f (Ny ) − Ny f Ny
N
2N 2
y∈SY
y∈SY


 ′ 2 
′
1
E Ny f (Ny ) − Ny f Ny
= O
N2


1
V [Ny f (Ny )] ,
= O
N2

(32)

′

where the second step follows from the fact that SY is finite and the last step follows from the fact that Ny and Ny are iid.
The expression V [Ny f (Ny )] is simply a sum of terms of the form of V Nγy where 0 < γ ≤ 1. Even the covariance terms
can be bounded by the square root of the product of these terms by the Cauchy Schwarz inequality.
Let py = fY (y). Consider the Taylor series expansion of the function h(x) = xγ at the point N py . This is
h(x)

=

γ

γ−1

(N py ) + γ (N py )

(x − N py ) +

γ(γ − 1)
γ−2
2
(N py )
(x − N py )
2

∞
X
γ(γ − 1) . . . (γ − k + 1)
γ−k
2
(N py )
(x − N py ) .
(33)
k!
k=3

From Riordan
[36], we know that the ith central moment of Ny is O N ⌊i/2⌋ . Then since γ ≤ 1, the last terms in (33) are

O N −1 when x = Ny and we take the expectation. Thus

+

 
E Nγy

 2
=⇒ E Nγy


γ(γ − 1)
γ−1
(N py )
(1 − py ) + O N −1
2
2

γ(γ − 1)
2γ−2
2γ
2γ−1
(N py )
= (N py ) + γ(γ − 1)(1 − py ) (N py )
+
2

+O N −1 .
γ

= (N py ) +

By a similar Taylor series expansion, we have that



2γ
2γ−1
= (N py ) + γ(2γ − 1)(1 − py ) (N py )
+ O N −1 .
E N2γ
y
Combining these results gives

 
V Nγy


 2

− E Nγy
= E N2γ
y

= O N 2γ−1 + N 2γ−2 + N −1
= O (N ) ,



where the last step follows from the fact that γ ≤ 1. Combining this result with (32) gives
 
ii
h h
1
V E G̃hX ,hX|Y Y1 , . . . , YN = O
.
N
i
h

By the law of total variance, V G̃hX ,hX|Y = O N1 .
C. Extension to the Generalized Case

In this section, we sketch the theory for the case where both X and Y have a mixture of discrete and continuous components.
Denote the discrete and continuous components of X as X1 and X2 , respectively. Similarly, denote the discrete and continuous
components of Y as Y1 and Y2 , respectively. Then the generalized mutual information is
X Z  fX (x1 , x2 )fY (y1 , y2 ) 
fXY (x1 , x2 , y1 , y2 )dx2 dy2
G(X; Y) =
g
fXY (x1 , x2 , y1 , y2 )
y1 ∈SY1
x1 ∈SX1

=

X

y1 ∈SY1
x1 ∈SX1

fX1 Y1 (x1 , y1 )

Z

g



fX1 (x1 )fY1 (y1 )fX2 |X1 (x2 |x1 )fY2 |Y1 (y2 |y1 )
fX1 Y1 (x1 , y1 )fX2 Y2 |X1 Y1 (x2 , y2 |x1 y1 )



fX2 Y2 |X1 Y1 (x2 , y2 |x1 y1 )dx2 dy2 .

P
Define Ny1 = N
i=1 1{Y1,i =y1 } where Y1,i is the discrete component of Yi . Then the estimator we use for fY1 (y1 ) is Ny1 /N .
The estimators for fX1 (x1 ) and fX1 Y1 (x1 , y1 ) are defined similarly with Nx1 and Nz1 .

We first consider the conditional bias of the resulting plug-in estimator where we condition on the discrete random variables.
Recall that by Taylor series expansions, we decompose the bias into “variance-like” terms in (23) and “bias-like” terms in (19).
For the bias-like term, if we condition on the discrete random variables, then
in (20) in this case is

 the equivalent expression
i

N

multiplied by Nx1 /N . This results in terms of the form of, for example, Nz1 − fZ1 (x1 , y1 ) . The expected value of these
expressions is the ith central moment of a binomial random variable divided by N i which is O(1/N ) for i ≥ 1. Thus these

i
N
terms contribute O(1/N ) to the bias. In all other cases, the expected value of Nz1
is O(1). Thus only the constants are
affected by these terms in the equivalent expression in (22). Similar results hold for the estimators of fX1 (x1 ) and fY1 (y1 ).
For the “variance-like” terms, we can simply factor out the estimators for fX1 (x1 ), fY1 (y1 ), and fZ1 (z1 ). The expected
value of these estimators is again O(1) so they only affect the constants.
For the variance, the law of total variance can again be used by conditioning on the discrete components. For the conditional
variance, the Lipschitz conditions on g in this case simply scales the resulting terms by the square of the estimators for
fX1 (x1 ), fY1 (y1 ), and fZ1 (z1 ). Then since the expected value of the square of these estimators is O(1), the expected value of
the conditional variance is still O(1/N ). Then by similar arguments given above for the bias and in Section E-B, the variance
of the conditional expectation of the estimator is also O(1/N ). Thus the total variance is O(1/N ).
P ROOF

A PPENDIX F
T HEOREM 6 (CLT)

OF

This proof shares some similarities with the CLT proof for the divergence functional estimators in [27], [28]. The primary
differences again deal with handling products of marginal density estimators and with handling two of the terms in the
Efron-Stein inequality. We will first find the asymptotic distribution of
"
!
!#!
N
i
h
√ 
f̃X,hX (Xi )f̃Y,hY (Yi )
1 X
f̃X,hX (Xi )f̃Y,hY (Yi )
g
=√
N G̃hX ,hY − E G̃hX ,hY
− EZi g
N i=1
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )
"
!#
"
!#!
N
f̃X,hX (Xi )f̃Y,hY (Yi )
f̃X,hX (Xi )f̃Y,hY (Yi )
1 X
−E g
.
EZi g
+√
N i=1
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )
By the standard central limit theorem [40], the second term converges in distribution to a Gaussian random variable with
variance
!##
" "
f̃X,hX (X)f̃Y,hY (Y)
.
V EZ g
f̃Z,hZ (X, Y)
All that remains is to show that the first term converges in probability to zero as Slutsky’s theorem [41] can then be applied.
Denote this first term as WN and note that E [WN ] = 0.
We will use Chebyshev’s n
inequality combined
o with the Efron-Stein inequality to bound ′ the variance of WN . Consider the
′
samples {Z1 , . . . , ZN } and Z1 , Z2 , . . . , ZN and the respective sequences WN and WN . This gives
"
!
!#!
′
1
f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃X,hX (X1 )f̃Y,hY (Y1 )
WN − WN = √
− EZ1 g
g
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )
N
!
"
!#!
′
′
′
′
f̃X,hX (X1 )f̃Y,hY (Y1 )
1
f̃X,hX (X1 )f̃Y,hY (Y1 )
− EZ′ g
+√
g
′
′
′
′
1
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )
N
!
!!
′
′
N
f̃X,hX (Xi )f̃Y,hY (Yi )
1 X
f̃X,hX (Xi )f̃Y,hY (Yi )
+√
−g
.
(34)
g
′
N i=2
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )
Note that

E g

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

!

"

− EZ1 g

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

!#!2 

"

"

 = E VX1 g

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

!##

.

h 
i
f̃
(X )f̃
(Y1 )
Y
We will use the Efron-Stein inequality to bound VX1 g X,hf̃X 1(XY,h
. We thus need to bound the conditional
1 ,Y1 )
Z,hZ
expectation of the term
!
!2
′
′
f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃X,hX (X1 )f̃Y,hY (Y1 )
−g
,
g
′
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )

′

where Zi is replaced with Zi in the KDEs for some i 6= 1. Using similar steps as in Section D, we have that

!
! 2
′
′


f̃
(X
)
f̃
(Y
)
1
(Y
)
(X
)
f̃
f̃
1
1
1
1
Y,h
X,h
X,h
Y,h
Y
X
X
Y


.
−g
=O
E g
′
N2
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )

h 
i
f̃
(X )f̃
(Y1 )
Y
Then by the Efron-Stein inequality, VX1 g X,hf̃X 1(XY,h
=O
,Y )
Z,hZ



E

1
N

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

g

A similar result holds for the g
For the third term in (34),


E



′

!

′

f̃X,hX (X1 )f̃Y,hY (Y1 )
′

′

f̃Z,hZ (X1 ,Y1 )

N
X

g

i=2

1

1

"


. Therefore

f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃Z,hZ (X1 , Y1 )

− EZ1 g


1
N

!#!2 

=O



1
N2



.

term in (34).

f̃X,hX (Xi )f̃Y,hY (Yi )
f̃Z,hZ (Xi , Yi )

!

′

′

f̃X,hX (Xi )f̃Y,hY (Yi )

−g

′

f̃Z,hZ (Xi , Yi )

! !2 


!
!
′
′
f̃X,hX (Xi )f̃Y,hY (Yi )
f̃X,hX (Xi )f̃Y,hY (Yi )
−g
E g
=
′
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )
i,j=2
!
!#
′
′
f̃X,hX (Xj )f̃Y,hY (Yj )
f̃X,hX (Xj )f̃Y,hY (Yj )
−g
× g
′
f̃Z,hZ (Xj , Yj )
f̃Z,hZ (Xj , Yj )
N
X

"

For the N − 1 terms where i = j, we know from Section D that

!
! 2
′
′


f̃
(X
)
f̃
(Y
)
(Y
)
(X
)
f̃
f̃
1
i
i
i
i Y,hY
X,hX
X,hX
Y,hY

.
E g
−g
=
O
′
N2
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )

Thus these terms contribute O(1/N ). For the N 2 − N terms where i 6= j, we can do multiple substitutions of the form
X −X
uj = jhX 1 resulting in
"
!
!
′
′
f̃X,hX (Xi )f̃Y,hY (Yi )
f̃X,hX (Xi )f̃Y,hY (Yi )
E g
−g
′
f̃Z,hZ (Xi , Yi )
f̃Z,hZ (Xi , Yi )
!
!#
!
′
′
X 2dY
f̃X,hX (Xj )f̃Y,hY (Yj )
h2d
f̃X,hX (Xj )f̃Y,hY (Yj )
X hY
−g
× g
=O
.
′
N2
f̃Z,hZ (Xj , Yj )
f̃Z,hZ (Xj , Yj )
Since hdXX hdYY = o(1),


E

N
X
i=2

g

f̃X,hX (Xi )f̃Y,hY (Yi )
f̃Z,hZ (Xi , Yi )

!

′

−g

′

f̃X,hX (Xi )f̃Y,hY (Yi )
′

f̃Z,hZ (Xi , Yi )

! !2 

 = o(1).

Combining all of these results with Jensen’s inequality gives

"
!
!#!2 



2
′
f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃X,hX (X1 )f̃Y,hY (Y1 )
3

− EZ1 g
E WN − WN
≤ E g
N
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )

!
"
!#!2 
′
′
′
′
f̃X,hX (X1 )f̃Y,hY (Y1 )
f̃X,hX (X1 )f̃Y,hY (Y1 )
3 

− EZ′ g
g
+ E
′
′
′
′
1
N
f̃Z,hZ (X1 , Y1 )
f̃Z,hZ (X1 , Y1 )

!
!!!2 
′
′
N
f̃X,hX (Xi )f̃Y,hY (Yi )
3  X
f̃X,hX (Xi )f̃Y,hY (Yi )

+ E
−g
g
′
N
(X
,
Y
)
(X
,
Y
)
f̃
f̃
i
i
i
i
Z,h
Z
Z,hZ
i=2
 
1
.
=o
N

Applying the Efron-Stein inequality gives that V [WN ] = o(1). Then by ChebyShev’s inequality, WN converges to zero in
probability. This completes the proof for the plug-in estimator.
For the weighted ensemble estimator, we can write
!
N
i
h
X
√ 
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
1 X
N G̃w − E G̃w = √
w(lX , lY ) g
N i=1 lX ∈LX ,lY ∈LY
f̃Z,hZ (l) (Xi , Yi )
"
!#!
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
−EZi g
f̃Z,hZ (l) (Xi , Yi )


!
N
X
X
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
1

EZi 
w(lX , lY )g
+√
f̃Z,hZ (l) (Xi , Yi )
N i=1
lX ∈LX ,lY ∈LY

!
X
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
 .
−E 
w(lX , lY )g
f̃Z,hZ (l) (Xi , Yi )
lX ∈LX ,lY ∈LY

By the central limit theorem, the second term converges in distribution to a zero-mean Gaussian random variable with variance


!
X
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
 .
w(lX , lY )g
V EZi 
f̃
(X
,
Y
)
i
i
Z,h
(l)
Z
lX ∈LX ,lY ∈LY

From the previous results, the first term converges to zero in probability as it can be written as
!
N
X
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
1 X
w(lX , lY ) √
g
f̃Z,hZ (l) (Xi , Yi )
N i=1
lX ∈LX ,lY ∈LY
"
!#!
X
f̃X,hX (l) (Xi )f̃Y,hY (l) (Yi )
w(lX , lY )oP (1)
−EZi g
=
f̃Z,hZ (l) (Xi , Yi )
lX ∈LX ,lY ∈LY
= oP (1),

where oP (1) denotes convergence to zero in probability and we use the fact that linear combinations of random variables that
converge in probability individually to constants converge in probability to the linear combination of the constants. The proof
is finished with Slutsky’s theorem.
Note that the proof of Corollary 7 follows a similar procedure as the extension to the ensemble case.

