Ensemble Estimation of Generalized Mutual Information with Applications to Genomics
Kevin R. Moon Kumar Sricharan Alfred O. Hero III

#1𝐞̃_#1,h_#1
𝐆̃_h_X,h_X|y
𝐆̃_h_X(l),h_Y(l)
𝐆̃_h_X,h_Y
𝐆̃_h
𝐆̂_k_1,k_2
𝐆̂_k
𝐆̂_k(l)
#1𝐆̃_#1
Cov

#1𝐟̃_#1,h_#1
#1𝐟̃_#1,h
#1𝐟̃_#1,h_#1(l)
#1𝐞̃_#1,h_#1
ƪ#1𝐞̃_#1,h
#1#2𝐟̃_#1,h_#1(l_#2)

𝔼
𝕍
𝔹
𝒮
#1𝔼_#1

#1ρ_#1,k_#1
#1ρ_#1,k_#1+1

#1𝐟̂_#1,k_#1
#1#2𝐟̂_#1,k(#2)
#1𝐟̂_#1,k_#1+1
#1𝐟̅_#1,k_#1
#1𝐟̅_#1,k_#1+1

#1#2𝐞̃_#1#2,h_#1,h_#2
𝐘
𝐗
𝐙
𝐍
𝐖
𝐡_X_C|x
𝐡_Y_C|y

#1#2#3#1_#2^(#3)








Ensemble Estimation of Generalized Mutual Information with Applications
to Genomics
    Kevin R. Moon1, Kumar Sricharan2, Alfred O. Hero III3

1Dept. of Mathematics and Statistics, Utah State University,  kevin.moon@usu.edu 

2Intuit Inc., sricharankumar@intuit.com 

3EECS Dept., University of Michigan, hero@eecs.umich.edu  This work was supported in part by the US Army Research Office under grants W911NF1910269 and W911NF1510479, and by the National Nuclear Security Administration in the US Department of Energy under grant DE-NA0003921. This paper appeared in part in the Proceedings of the 2017 IEEE Intl. Symposium on Information Theory (ISIT) <cit.>.
    December 30, 2023
================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


Mutual information is a measure of the dependence between random variables
that has been used successfully in myriad applications in many fields.
Generalized mutual information measures that go beyond classical Shannon
mutual information have also received much interest in these applications.
We derive the mean squared error convergence rates of kernel density-based
plug-in estimators of general mutual information measures between
two multidimensional random variables  and  for two cases:
1)  and  are continuous; 2)  and  may have a mixture
of discrete and continuous components. Using the derived rates, we
propose an ensemble estimator of these information measures called
GENIE by taking a weighted sum of the plug-in estimators with varied
bandwidths. The resulting ensemble estimators achieve the 1/N parametric
mean squared error 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 the
mixture case, which frequently arises in applications (e.g. variable
selection in classification). The estimator is simple to implement
and 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 estimators and minimax rates are derived for the
continuous case. We demonstrate the ensemble estimator for the mixed
case on simulated data and apply the proposed estimator to analyze
gene relationships in single cell data.



mutual information; nonparametric estimation; central limit theorem; single cell data; feature selection; minimax rate




§ INTRODUCTION


Mutual information (MI) is a measure of the amount
of shared information between a pair of random variables  and
. MI estimation is related to the problem of estimating functionals
of probability distributions, which has received deserved attention
in recent years <cit.>.
Many statistical problems rely in some form upon accurate estimation
of functionals of probability distributions including estimating the
decay rates of error probabilities <cit.>, estimating
bounds on the Bayes error rate <cit.>,
and hypothesis testing <cit.>.
MI estimation, in particular, also has many applications in information
theory and machine learning including independent subspace analysis <cit.>,
structure learning <cit.>, fMRI data processing <cit.>,
forest density estimation <cit.>, clustering <cit.>,
neuron classification <cit.>, blind source
separation <cit.>, intrinsically motivated reinforcement
learning <cit.>, as well
as other data science applications such as sociology <cit.>,
computational biology <cit.>,
and improving neural network models <cit.>. A particularly
common application is feature selection or extraction where features
are chosen to maximize the MI between the chosen features (represented
by 𝐗) and the outcome variables (represented by 𝐘)
<cit.>. 

In many of these applications, the variables  and  may have
any mixture of discrete and continuous components. In feature selection,
for example, the predictor labels may have discrete components (e.g.
classification labels) while the input variables may have a mixture
of discrete and continuous features. 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 (N is the number of samples) in this setting where 
and/or  contain a mixture of discrete and continuous components.
Instead, most existing estimators of MI focus on the cases where both
 and  are either purely discrete or purely continuous. Also,
while many nonparametric estimators of MI exist, most have not been
generalized beyond Shannon or Rényi information. Furthermore,
minimax convergence rates are currently unknown for the continuous
and the mixture cases.

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. This framework can be applied
to accurately estimate general MI measures when either  and 
are purely continuous or the mixed case when  and  may contain
a mixture of discrete and continuous components. We derive an MI estimator
for these cases that achieves the parametric MSE rate when the conditional
densities of the continuous variables are sufficiently smooth, thus
achieving the minimax rate (which we also derive) in this setting.
We call this estimator the Generalized ENsemble
Information Estimator (GENIE).

Our estimation method applies to other MI measures in addition to
Shannon information, which have been the focus of much recent interest.
An information measure based on a quadratic divergence was defined
in <cit.>. A density-resampled version of MI
was introduced in <cit.> to better measure
gene relationships in single-cell data when sampling may not be uniform.
A MI measure based on the Pearson divergence was considered in <cit.>.
Minimal spanning tree <cit.> and generalized nearest-neighbor
graph <cit.> approaches have been developed for
estimating Rényi information <cit.>,
which has been used in many applications (e.g. <cit.>).



 §.§ Related Work


Many estimators for MI have been previously developed. Nearly all
of these estimators ignore the mixed case and focus on the case where
both  and  are either purely continuous or purely discrete.
A popular k-nearest neighbor (nn)-based estimator was proposed
in <cit.> which is a modification of the entropy
estimator derived in <cit.>. However, these
estimators have only been shown to achieve the parametric convergence
rate when the dimension of each of the random variables is less than
3 <cit.>. Furthermore, these estimators focus
only on estimating the Shannon MI between purely continuous random
variables. Similarly, the estimators in <cit.>
do not achieve the parametric rate and focus on the purely continuous
case. An adaptation of the Shannon MI estimator in <cit.>
was recently proposed to handle the discrete-continuous mixture case <cit.>.
While this estimator has been proven to be consistent, its convergence
rate is currently unknown. Central limit theorems have also been derived for several entropy estimators <cit.>, which can then be applied to Shannon MI. However, it is not clear if these results can be extended to more general MI functionals.

A neural network-based estimator of Shannon MI was proposed in <cit.>.
While this estimator is computationally efficient, its statistical
properties are largely unknown as the authors only prove convergence
in probability rates. It is also unclear how to extend this estimator
to other MI measures such as the Rényi information. A jackknife
approach to estimating Shannon MI was also recently proposed <cit.>.
This approach provides an automatic selection of the kernel bandwidth
for a plug-in kernel density estimator (KDE) and does not require
boundary correction, which is generally a major difficulty in estimating
functionals of probability distributions. However, the MSE convergence
rate of this estimator is also unknown. 

Much work has focused on the problem of estimating the entropy of
purely discrete random variables <cit.>.
Shannon MI can then be estimated by estimating the joint and marginal
entropies of  and . However, it is not clear if discrete
methods can be extended successfully to the mixed-case. Quantizing
the continuous components of the data is one potential approach that
has been shown to be consistent for some quantization schemes in the
purely continuous case <cit.> but it is currently
unknown if similar approaches can be applied in the mixed-case. Also,
extending these estimators to general MI measures like Rényi
information is not straightforward.

Recent work has focused on nonparametric divergence estimation for
continuous random variables. One approach <cit.>
uses an optimal KDE to achieve the parametric convergence rate when
the densities are at least d <cit.>
or d/2 <cit.>
times differentiable where d is the dimension of the data. These
methods, like ours, assume that the densities are bounded away from
zero as this simplifies the analysis. However, this induces a boundary
on the densities' support set. For accurate estimation, the optimal
KDE approaches 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.

Some methods for estimating distributional functionals have relaxed
the boundedness assumption on the densities <cit.>.
These approaches typically assume that the tails of the densities
decay at a sufficiently fast rate (e.g. sub-exponential or sub-Gaussian).
In <cit.>, the authors only consider densities
with up to 2 derivatives as it is  difficult to exploit higher smoothness when the densities are not lower-bounded.

More closely related work <cit.>
uses an ensemble approach to estimate entropy or divergence functionals
for continuous random variables. These works construct an ensemble
of simple plug-in estimators by varying the neighborhood size of 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 <cit.>
or d/2 <cit.> 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. The ensemble approach also automatically corrects
for bias at the boundary of the densities' support set.

Finally, <cit.> 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 this isn't an issue when the densities
are smooth <cit.>.

For the mixture setting, we focus on the important special case where
the components of each observation are assumed to decompose into discrete
and continuous dimensions. This enables the density to be factored:
f_X(x)=f_X_D(x_D)f_X_C|X_D(x_C|x_D)
where x_C and x_D are the continuous and discrete components
of x. We note that this excludes the more general case considered
by <cit.> where one or more components can have
discrete and continuous values simultaneously. However, our setting
is a common occurrence in many machine learning and statistical problems.
For example, a search within the UCI Machine Learning Repository <cit.>
yields many datasets with such structure. Many statistical models
have also focused on similar settings <cit.>.
Thus we believe that this special case warrants its own treatment
and retain the more general case for future work. Despite the importance
of this mixed setting, no other MI estimators have been derived or
analyzed that achieve the parametric MSE convergence rate.



 §.§ Contributions


In the context of this related work, we make the following novel contributions
in this paper: 


  * For purely continuous random variables, we derive the asymptotic bias
and variance of kernel density plug-in MI estimators for general MI
measures without boundary correction <cit.>
(Section <ref>). 

  * We leverage the results for the purely continuous case to derive the
bias and variance of general kernel density plug-in MI estimators
when  and/or  contain a mixture of discrete and continuous
components by reformulating the densities as a mixture of the conditional
density of the continuous variables given the discrete variables (Section
<ref>). Note that this is a special case of the mixture
setting where discrete and continuous components are separated into
different dimensions.

  * We leverage this theory for the mixed case described above in conjunction
with the generalized theory of ensemble estimators <cit.>
to derive GENIE. To the best of our knowledge, this is the first non-parametric
estimator of general MI measures that achieves a parametric rate of
MSE convergence of O(1/N) when the densities are sufficiently smooth  for any mixed case (Section
<ref>) let alone the special case we consider,
where N is the number of samples available from each distribution.

  * We prove a minimax lower bound for the convergence rate of MI estimators
in the purely continuous case (Section <ref>). This
unifies the minimax theory for estimating continuous entropy <cit.>
and divergence functionals <cit.>.
Neither of these approaches are directly extendable to the MI case
due to the dependence of the marginal distributions on the joint distribution
and the integral relationship between the joint and the marginals.
Therefore, we have tailored the proof to the MI estimation case. We
also show that the MI ensemble estimator achieves the minimax rate
when the densities are sufficiently smooth.

  * We derive a central limit theorem for the ensemble estimators (Section
<ref>).

  * We apply the method to single-cell RNA-sequencing feature selection
problems (Section <ref>).

We note that KDE plug-in approaches to estimating functionals such
as entropy and MI are well-known and perhaps the simplest approach <cit.>.
Applying the generalized theory of ensemble estimation to the KDE
plug-in estimator does not raise the complexity of the estimators
substantially, either computationally or conceptually. Yet by employing
these simple methods, the resulting ensemble estimator is able to
achieve the minimax convergence rate for sufficiently smooth densities
without employing more complicated von-Mises expansions (as in <cit.>)
or boundary correction (as in <cit.>)
to reduce the bias.



§ MUTUAL INFORMATION FUNCTIONALS


We first define a family of MI functionals based on f-divergence
functionals which are defined as follows. Let P and Q be probability
measures on the Euclidean space 𝒮. Let g:(0,∞)→ℝ be the f-divergence shaping function.
The f-divergence functional associated with g is <cit.>

    D_g(P||Q):=_Q[g(dP/dQ)],

where dP/dQ is the Radon-Nikodym derivative and _Q
indicates the expectation wrt to the measure Q. To obtain a true
divergence, we require g to be convex and g(1)=0. However, we
consider more general functionals and so we do not place these restrictions
on g.

A generalized MI functional can be derived from (<ref>).
Let  and  be (potentially multivariate) random variables
with respective marginal probability measures P_X and P_Y
and joint probability measure P_XY. Let g be as before. Then
the MI functional associated with g is 

    I(;):=D_g(.P_XP_Y‖ P_XY).

Shannon MI can be obtained from (<ref>) by setting
g(t)=-log t.

If  and  are purely continuous random variables with respective
marginal probability densities f_X and f_Y and joint probability
density f_XY, then (<ref>) can be written as 

    I(;)=∫ g(f_X(x)f_Y(y)/f_XY(x,y))f_XY(x,y)dxdy.


However, we are also interested in the case where  or  may
have a mixture of discrete and continuous components. In this special
case, the distributions can be factored into a product of the conditional
density and the probability mass functions. The MI can then be expressed
as a sum of integrals which can then be individually estimated. To
do this, denote the continuous and discrete components of  as
_C and _D, respectively. Denote _C and _D
similarly. Let =(𝐗,𝐘)^T and let
_C and _D be the respective continuous and discrete
components of . Consider the probability distributions f_XY,
f_X, f_Y and the corresponding densities that are obtained
by conditioning on _D and _D, e.g. f_XY(x_C,x_D,y_C,y_D)=f_X_CY_C|X_DY_D(x_C,y_C|x_D,y_D)f_X_DY_D(x_D,y_D).
Then after factoring the distributions, (<ref>) can
be written as

    I(;)     =∑_x_D,y_D∫ g(f_X(x_C,x_D)f_Y(y_C,y_D)/f_XY(x_C,x_D,y_C,y_D))dF_XY(x_C,x_D,y_C,y_D)
        =∑_z_Df_Z_D(z_D)∫ g(R_1(z_C)R_2(z_D))f_Z_C|Z_D(z_C|z_D)dz_C,

where 

    R_1(z_C)     =f_X_C|X_D(x_C|x_D)f_Y_C|Y_D(y_C|y_D)/f_X_CY_C|X_DY_D(x_C,y_C|x_D,y_D),
    
    R_2(z_D)     =f_X_D(x_D)f_Y_D(y_D)/f_X_DY_D(x_D,y_D).

The expression R_1 is the ratio of the product of the conditional
densities f_X_C|X_D and f_Y_C|Y_D to the conditional
density f_X_CY_C|X_DY_D. It is a continuous function
of z_C. Similarly, the expression R_2 is the ratio of the
product of the probability mass functions (pmf) f_X_D and f_Y_D
to the pmf f_X_DY_D and is a discrete function of z_D. 

In the following sections, we will obtain MSE convergence rates of
KDE plug-in estimators of the general MI measures described above.
We first focus on the case when  and  are purely continuous
(Equation (<ref>)). We then generalize to the case where
 and  may have any mixture of continuous and discrete components
(Equation (<ref>)). The derived convergence rates can
then be used to derive ensemble estimators that achieve the parametric
MSE rate.



§ CONTINUOUS RANDOM VARIABLES


For this section, we define KDE plug-in estimators
of general MI measures under the assumption that  and  are
purely continuous. Thus _C= and _C= and we can write 


    I(;)=∫ g(f_X(x)f_Y(y)/f_XY(x,y))f_XY(x,y)dxdy.

We then derive the MSE convergence rate of the KDE plug-in estimator.
We also present a minimax lower bound for MI estimation in this continuous
setting.

To more easily generalize our results to the mixture case, we consider
a modified version of (<ref>) where the densities are
weighted as follows. Let ν be a 3-dimensional vector with 0<ν_i≤1
for each i∈{1,2,3}. We can then write 

    I_ν(;)=∫ g(f_X(x)f_Y(y)ν_1ν_2/f_XY(x,y)ν_3)f_XY(x,y)dxdy.

The expression in (<ref>) reduces to that in (<ref>)
when ν_i=1 for each i∈{1,2,3}. When we generalize to
the mixture case, the pmf estimators will be substituted into ν.



 §.§ The KDE Plug-in Estimator


Let f_X(x), f_Y(y), and f_XY(x,y) be d_X, d_Y,
and d_X+d_Y=d-dimensional probability densities. Since we are
assuming for now that 𝐗 and 𝐘 are continuous
with marginal densities f_X and f_Y, the MI functional I_v(;)
can be estimated using KDEs. Assume that N i.i.d. samples {𝐙_1,…,𝐙_N}
are available from the joint density f_XY with 𝐙_i=(𝐗_i,𝐘_i)^T.
Let h_X, h_Y be kernel bandwidths. Let K_X(·) and
K_Y(·) be symmetric kernel functions with ∫ K_X(x)dx=∫ K_Y(y)dy=1,
||K_X||_∞, ||K_Y||_∞<∞ where ||K||_∞=sup_x|K(x)|.
The KDEs for f_X, f_Y, and f_XY=f_Z, respectively,
are 

    X(x)     =    1/Nh_X^d_X∑_i=1^NK_X(x-𝐗_i/h_X),
     Y(y)     =    1/Nh_Y^d_Y∑_i=1^NK_Y(y-𝐘_i/h_Y),
     Z(x,y)     =    1/Nh_X^d_Xh_Y^d_Y∑_i=1^NK_X(x-𝐗_i/h_X)K_Y(y-𝐘_i/h_Y),

where h_Z=(h_X,h_Y). Then I_ν(;) can be estimated
with a KDE plug-in estimator: 

    =1/N∑_i=1^Ng( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3).

Note that in this estimator we evaluate the KDEs at each of the data
points. In practice, this is done using a leave-one-out KDE. This
enables us to avoid evaluating a high-dimensional integral and instead
estimate the integral with the empirical average in eq. (<ref>).



 §.§ MSE Convergence Rate of the Continuous Plug-in Estimator


We are interested in the MSE convergence rate of the KDE plug-in estimator
in eq. (<ref>). The MSE of an estimator can be expressed
as the sum of the squared bias and the variance of the estimator.
We first focus on the bias of the estimator . The bias of nonparametric
estimators typically depends on the smoothness of the functions that
are being estimated. In our case, we have multiple functions including
the joint and marginal densities and the function g. We quantify
the smoothness of the densities using the Hölder class Σ(s,H):

 Let
𝒳⊂ℝ^d be a compact space. For q=(q_1,…,q_d),
q_i∈ℕ, define |q|=∑_i=1^dq_i and D^q=∂^|q|/∂ x_1^q_1…∂ x_d^q_d.
The Hölder class Σ(s,H) of functions on L_2(𝒳)
consists of the functions f that satisfy 

    |D^qf(x)-D^qf(y)|≤ H‖ x-y‖ ^s-|q|,

for all x, y∈𝒳 and for all q s.t. |q|≤⌊ s⌋.


A key fact that comes from Definition <ref> is that if
a function f belongs to Σ(s,H), then it is r=⌊ s⌋
times differentiable. Given this definition, the full assumptions
we make to derive bias convergence rates are:


  * (𝒜.0): The kernels K_X and K_Y are symmetric
product kernels with bounded support. 

  * (𝒜.1): There exist constants ϵ_0, ϵ_∞
such that 0<ϵ_0≤ f_X(x)≤ϵ_∞<∞
∀ x∈𝒮_X, ϵ_0≤ f_Y(y)≤ϵ_∞
∀ y∈𝒮_Y, and ϵ_0≤ f_XY(x,y)≤ϵ_∞
∀(x,y)∈𝒮_X×𝒮_Y.

  * (𝒜.2): Each of the densities belong to Σ(s,H)
in the interior of their support sets with s≥2.

  * (𝒜.3): g(t_1/t_2) has an infinite number of mixed
derivatives wrt t_1 and t_2. 

  * (𝒜.4): |∂^k+lg(t_1/t_2)/∂ t_1^k∂ t_2^l|/(k!l!),
k,l=0,1,… are strictly upper bounded for ϵ_0≤ t_1,t_2≤ϵ_∞.

  * (𝒜.5): Let K be either K_X or K_Y, 𝒮
either 𝒮_X or 𝒮_Y, h either h_X
or h_Y, f either f_X or f_Y,  and d either d_X
or d_Y. Let q=(q_1,…,q_d) with q_i∈ℕ and |q|≤ r=⌊ s⌋. Then we assume for any positive integers t and l that 

    ∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮K^l(u)u^q D^q f(x) du)^tdx = v_t(h),

 where v_t(h) admits the expansion 

    v_t(h)=∑_i=1^r-|q|e_i,q,t,lh^i+o(h^r-|q|),

 for some constants e_i,q,t,l.














 


These assumptions can largely be summarized as follows: 1) f_X,
f_Y, f_XY, and g are smooth (𝒜.2-𝒜.4)
; 2) f_X and f_Y have bounded support sets 𝒮_X
and 𝒮_Y with respective dimensions d_X and d_Y
(𝒜.1); 3) f_X, f_Y, and f_XY are strictly
lower bounded on their support sets (𝒜.1); and 4) the
boundary of the support set is smooth (𝒜.5). More specifically,
assumption 𝒜.5 states that the support of the density
is smooth with respect to the kernel K in the sense that the expected value of a polynomial with coefficients consisting of the densities and their derivatives near the boundary is a smooth function of the bandwidth h. The inner integral in (<ref>) captures this
expectation while the outer integral averages this inner integral
over all points near the boundary of the support. The v_t(h)
term captures the fact that the smoothness of this expectation is
proportional to the smoothness of the function D^q f(x). 












While these assumptions may appear highly technical, they are satisfied
for relatively simple support sets and for common kernels, functions
g, and densities and thus are widely applicable <cit.>.
These assumptions are also comparable to those in similar studies
on asymptotic convergence analysis <cit.>.
Some studies consider the case where the densities are not strictly lower bounded, which makes the problem different <cit.> with different minimax rates (see <cit.> for the entropy estimation case).


In particular, assumption 𝒜.5  is satisfied if the kernel K is smooth, has either circular or rectangular support (which includes product kernels), and the density support set consists of the unit cube. See Appendix <ref> for details. The unit cube assumption is common in the nonparametric density functional estimation literature <cit.> as the results can then be extended to density support sets that are isomorphic to the unit cube. 

To derive the convergence rates of many state-of-the art distributional functional estimators, authors commonly assume that the derivatives of the density f(x) vanish near the boundary <cit.>. Note that in this assumption, the density itself is not required to vanish near the boundary. Thus densities such as the uniform distribution satisfy this common assumption. However, this assumption is stronger than 𝒜.5 as formalized in Proposition <ref> below. Our weaker assumption 𝒜.5 comes at a small cost as we require the f-divergence shaping function g to be infinitely differentiable. In contrast, the authors in <cit.> assume that the shaping function has a finite number of derivatives. In practice, this tradeoff does not have a major practical impact as most shaping functions of interest are either infinitely differentiable everywhere (e.g. Shannon and Renyi information) or not differentiable on a set of measure zero (e.g. the total variation distance and the Bayes error rate in the divergence case).


Let the density support set 𝒮 be the unit cube. Let the derivatives of the density f up to order r vanish at the boundary of the density support set. Assume that ||K||_∞<∞ and the support of K is bounded with either rectangular or circular support. Then assumption 𝒜.5 is satisfied.


The proof is given in Appendix <ref>. The assumption of vanishing density derivatives at the boundary is strictly weaker than assumption 𝒜.5. As an example, consider a standard Gaussian distribution truncated to the support [-1,1]^d. Clearly, the derivatives of this density do not vanish at the boundary. However, we show in Appendix <ref> that this density satisfies 𝒜.5.

We note that the boundary assumption 𝒜.5 does not directly
result in parametric convergence rates for the plug-in estimator ,
which is in contrast with the boundary assumptions in <cit.>.
The estimators in <cit.>
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  without boundary correction, greatly simplifying
our estimator.


 Under assumptions 𝒜.0-𝒜.5,
the bias of  is 

    []     =    ∑_j=0
    
    i+j≠0
    ^r∑_i=0^rc_10,i,j(ν_1ν_2,ν_3)h_X^ih_Y^j+c_11/Nh_X^d_Xh_Y^d_Y
           +O(h_X^s+h_Y^s+1/Nh_X^d_Xh_Y^d_Y),

where the constants in (<ref>) are independent of the bandwidths
h_X and h_Y and depend on the
densities and their derivatives, the functional g and its derivatives,
and the kernels. They also include polynomial terms of ν_1ν_2
and ν_3 when ν_i≠ 1.
 

Expressions for
the constants in  (<ref>) are not given in this paper due to their complexity. These constants are not needed as the bias rates in Theorem <ref> are sufficient to implement ensemble bias reduction. The resultant ensemble estimator  achieves the parametric MSE convergence
rate O(1/N) (see Section <ref> for the mixed case and
Appendix <ref> for the continuous case). 

We also derive a refined expression for the bias that enables
us to achieve the parametric convergence rate under less restrictive
smoothness assumptions on the densities (s>(d_X+d_Y)/2 compared
to s≥ d_X+d_Y for (<ref>)). However, the resulting
expansion has more terms and the ensemble estimator is more complicated to implement. Thus we have chosen
to present the simpler case here. The more complex expansion and estimator are presented
in Appendix <ref>.

Having obtained an expression for the bias of , we now present
an upper bound on its variance to complete the derivation of its MSE.
 If the functional g
is Lipschitz continuous in both of its arguments with Lipschitz constant
C_g, then the variance of  is 

    []≤22C_g^2||K_X· K_Y||_∞^2/N.


The Lipschitz assumption on g for the variance result is comparable
to assumptions made by others for nonparametric estimation of distributional
functionals <cit.>
and is satisfied for Shannon and Renyi informations when the densities
are bounded above and below. Note that Theorem <ref>
requires much less strict assumptions than Theorem <ref>.
The proofs of Theorems <ref> and <ref> are
given in Appendix <ref> and <ref>,
respectively.

Theorems <ref> and <ref> indicate that for
the MSE to go to zero, we require h_X,h_Y→0 and Nh_X^d_Xh_Y^d_Y→∞.
In Section <ref>, we will use Theorems <ref>
and <ref> to derive bias and variance expressions for
the MI plug-in estimators under the more general cases where 
and/or  may contain a mixture of discrete and continuous components.
We will then use these convergence rate results to derive MI ensemble
estimators for both cases (purely continuous random variables and
mixed random variables) that achieve the parametric MSE convergence
rate regardless of the dimension as long as the densities are sufficiently
smooth. 



 §.§ Minimax Rate for MI estimation


We wrap up this section with a minimax lower
bound on the MSE rate of convergence for the continuous MI estimation
problem. 
 Assume that g is
at least twice differentiable and that given ϵ>0, |g”(ϵ)|>0.
Define the set of functions Σ(s,H,ϵ_0,ϵ_∞)
to be the set of Hölder continuous functions Σ(s,H) that
are bounded between ϵ_0 and ϵ_∞. Then
with γ=min{ 8s/(4s+d_X+d_Y),1},
there exists a strictly positive constant c such that 

    lim inf_N→∞inf_Ĝ_Nsup_f_XY∈Σ(s,H,ϵ_0,ϵ_∞)[(Ĝ_N-I(;))^2]≥ cN^-γ.


The proof uses Le Cam's method <cit.> and is
given in Appendix <ref>. Theorem <ref>
indicates that the minimax rate is the parametric rate N^-1 as
long as s≥(d_X+d_Y)/4. This is consistent with
minimax rates for divergence <cit.>
and entropy <cit.> functional estimation, thus
expanding the previous theory on minimax estimation of information
theoretic functionals. 

In Section <ref> and Appendix <ref>,
we derive MI estimators that achieve the minimax rate when s≥ d_X+d_Y and s>(d_X+d_Y)/2, respectively.
While estimators have been derived for the divergence estimation problem
that achieve the minimax rate for less smooth densities, they require
numerical integration and are thus computationally slow <cit.>.
Deriving estimators of these functionals (e.g. MI and divergence)
that are known to achieve the minimax rate in this less smooth regime
and that are computationally reasonable thus remains an open problem.






§ MIXED RANDOM VARIABLES


In this section, we extend the results of Section <ref>
to general MI estimation when 𝐗 and 𝐘 may
have a mixture of discrete and continuous components. We focus on
the most complex case:  and  both have discrete and continuous
components. The MI between  and  is written in (<ref>).



 §.§ KDE Plug-in Estimator


We first define the KDE plug-in estimator of (<ref>).
Let 𝒮_Y_C and 𝒮_X_C be the respective
supports of the corresponding densities of _C and _C
and let 𝒮_Y_D and 𝒮_X_D be the respective
supports of the corresponding probability mass functions of _D
and _D. Suppose we have N i.i.d. samples of (,) drawn
from f_XY where the ith samples are denoted as (_i,_i)=(_i,C,_i,D,_i,C,_i,D).
Define the following random variables:

    𝐍_y    =∑_i=1^N1_{𝐘_i,D=y} ,
    _x    =∑_i=1^N1_{_i,D=x} ,
    _xy    =∑_i=1^N1_{_i,D=x,_i,D=y} ,

where x∈𝒮_X_D, y∈𝒮_Y_D, and
1_{·} is the indicator function. These will be used to
estimate the pmfs of the discrete components of  and .

For the continuous components, we will condition on the discrete components
and construct KDEs for the conditional probability density functions.
Let 𝒮_X_C and 𝒮_Y_C be the respective
supports of the marginal densities f_X_C and f_Y_C with
corresponding dimensions of d_X and d_Y. As before, let
K_X(·) and K_Y(·) be kernel functions with ∫ K_X(x)dx=∫ K_Y(y)dy=1,
||K_X||_∞, ||K_Y||_∞<∞ where ||K||_∞=sup_x|K(x)|.
Consider the following sets:

    𝒳_x    ={._i,C∈{_1,C,…,_N,C}|_i,D=x} ,
    𝒴_y    ={._i,C∈{_1,C,…,_N,C}|_i,D=y} .

The set 𝒳_x is the set of the continuous  data
points where the corresponding discrete component is equal to x.
The set 𝒴_y is defined similarly. The KDEs for f_X_C|X_D,
f_Y_C|Y_D, and f_X_CY_C|X_DY_Dat x∈𝒮_X_D
and y∈𝒮_Y_D are, respectively, 

    X_C|x(x)     =1/_xh_X_C|x^d_X∑_[ _j,C∈𝒳_x;     i≠ j ]K_X(x-_j,C/h_X_C|x),
    Y_C|y(y)     =1/_yh_Y_C|y^d_Y∑_[ _j,C∈𝒴_y;     i≠ j ]K_Y(y-_j,C/h_Y_C|y),
    Z_C|z(x,y)     =1/_xyh_X_C|x^d_Xh_Y_C|y^d_Y∑_[ _j,C∈𝒴_y AND _j,C∈𝒳_x;                  i≠ j ]K_X(x-_j,C/h_X_C|x)K_Y(y-_j,C/h_Y_C|y),

where _C=(_C,_C) and h_Z_C|z=(h_X_C|x,h_Y_C|y).
Note that we allow the bandwidths to depend on the discrete components
of 𝐗 and 𝐘. The reason for this is that the
bandwidth is generally chosen as a function of the number of data
points, which will differ for these conditional distributions as the
discrete components of 𝐗 and 𝐘 differ.

The MI I(;) can then be estimated by plugging in the conditional
KDEs. Note that the MI in eq. (<ref>) is written as
a weighted sum of integral functionals. We therefore first define
an intermediate estimator of the integral functionals:

    h_X_C|x,h_Y_C|y=1/_xy∑__C∈𝒳_xAND_C∈𝒴_yg(X_C|x(_C)Y_C|y(_C)/Z_C|z(_C,_C)×_x_y/N_xy).

Again in practice, we evaluate the KDEs at each of the data points
using a leave-one-out KDE, enabling us to avoid evaluating a high-dimensional
integral. We then define a plug-in KDE estimator of I(;):

    h_X_C|X_D,h_Y_C|Y_D=∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/Nh_X_C|x,h_Y_C|y.


The quality of the conditional density estimates in terms of bias
and variance depends on the choice of bandwidths h_X_C|x and
h_Y_C|y. That is, for the KDE X_C|x to converge in
MSE, it is necessary that h_X_C|x→0 and _xh_X_C|x^d_X→0
as _x→∞ (a similar result holds for h_Y_C|y) <cit.>.
Furthermore, we will see when we derive the bias and variance of h_X_C|X_D,h_Y_C|Y_D
that these conditions are also necessary for h_X_C|X_D,h_Y_C|Y_D
to converge in MSE. Thus, when deriving the MSE convergence rate of
h_X_C|X_D,h_Y_C|Y_D, we will assume that h_X_C|x
is a function of _x and h_Y_C|y is a function of _y.



 §.§ MSE Convergence Rates of the Mixed Plug-in Estimator


Here we derive the MSE convergence rate
of a plug-in estimator of MI when the random variables have a mixture
of discrete and continuous components. We will need the following:

Let _y, _x, and _xy
be defined as in (<ref>). Assume that their corresponding
probability mass functions are bounded away from zero. If α∈ℝ\{0,1}
and λ+β+γ∈ℝ\{0,1}, then 

    [_xy^α]     =(Nf_X_DY_D(x,y))^α+O(N^α-1)
    [_xy^λ_x^β_y^γ]     =N^λ+β+γ(f_X_DY_D(x,y))^λ(f_X_D(x))^β(f_Y_D(y))^γ+O(N^λ+β+γ-1).


The proof is in Appendix <ref> and uses the
generalized binomial theorem, Taylor series expansions, and known
results about the central moments of binomial random variables <cit.>.
Lemma <ref> provides key results on moments
of products of the binomial random variables _xy, _x,
and _y. These results can be used to derive the bias and variance
of a plug-in estimator of MI with mixed components in (<ref>)
as long as the bias and variance of the corresponding plug-in estimator
for the continuous weighted case in (<ref>) is known.
This is demonstrated in the following theorems for the KDE plug-in
estimator h_X_C|X_D,h_Y_C|Y_D.
 Assume that assumptions 𝒜.0-𝒜.5
hold with respect to the functional g, the kernels K_X and
K_Y, and the densities f_X_C|X_D, f_Y_C|Y_D
and f_X_CY_C|X_DY_D. Assume that |𝒮_X_D|,|𝒮_Y_D|<∞.
Assume that =l_X𝐍_x^-β and =l_Y_y^-α
with 0<β<1/d_X, 0<α<1/d_Y, and
l_X,l_Y>0. Then the bias of h_X_C|X_D,h_Y_C|Y_D
is 

    [h_X_C|X_D,h_Y_C|Y_D]     =∑_i,j=0
    
    i+j≠0
    ^rc_13,i,jl_X^il_Y^jN^-iβ-jα+O(N^-sα+N^-sβ+N^β d_X+α d_Y-1).


The constants depend on the underlying densities, the chosen kernels,
the functional g, and the probability mass functions and are independent
of l_X, l_Y, and N. Furthermore, these rates are asymptotically
tight.
 Assume that =l_X𝐍_x^-β
and =l_Y_y^-α with 0<β<1/d_X,
0<α<1/d_Y, β d_X+α d_Y≤1, and
l_X,l_Y>0. Assume that |𝒮_X_D|,|𝒮_Y_D|<∞.
If the shaping function g is Lipschitz continuous in both of its arguments,
then the variance of h_X_C|X_D,h_Y_C|Y_D is O(1/N).

These theorems provide the necessary information for applying the
theory of optimally weighted ensemble estimation to obtain MI estimators
with improved rates (see Section <ref>). 



 §.§ Proof Sketches of Theorems <ref> and <ref>


For Theorem <ref>, the proof splits the bias term
into two terms by adding and subtracting g(𝒯(,)_x_y/N_xy)
for each pair (x,y) where 𝒯(,) is independent
of the data samples and is defined in Eq. (<ref>). It can
be shown that the newly added term has bias O(1/N). The other term
is handled by conditioning on the discrete components of the data
samples to obtain the conditional bias terms [.h_X_C|x,h_Y_C|y|_1,D,…,_N,D,_1,D,…,_N,D]
for each pair (x,y). Theorem <ref> can then be applied
to each of these terms to obtain expressions of the random variables
_x, _y, and _xy with terms of the form given
in Lemma <ref>. Lemma <ref>
can be applied to these terms to obtain the final result, where care
is taken to ensure that all relevant terms have been handled properly.
The full proof is given in Appendix <ref>.

To prove Theorem <ref>, we use the law of total variance
to split the variance into two terms: the expected value of the variance
conditioned on the discrete components of the data samples and the
variance of the conditional expectation. Theorem <ref>
is applied to the conditional variance term. For the conditional expectation
term, we use results obtained in the proof of Theorem <ref>
combined with the Efron-Stein inequality <cit.>
to obtain expressions of the random variables _x, _y,
and _xy. Lemma <ref> can be applied
again to these terms to obtain the final result. The full proof is
given in Appendix <ref>.



§ ENSEMBLE ESTIMATION OF GENERALIZED MI


 

If no bias correction is performed, then Theorems <ref>
and <ref> show that the optimal bias rate of the KDE
plug-in estimators  and h_X_C|X_D,h_Y_C|Y_D
is O(1/N^1/(d_X+d_Y+1)), which converges very
slowly to zero when either d_X or d_Y are not small. Thus
the standard KDE plug-in estimators will perform poorly in  higher-dimensional
settings. We use the theory of optimally weighted ensemble estimation
developed in <cit.> to improve this rate. For brevity,
we focus on the case where  and  both contain a mixture
of discrete and continuous components. The purely continuous case
is described in Appendix <ref>.

An ensemble of estimators is first formed by choosing different bandwidth
values for the plug-in estimators as follows. Let ℒ be
a set of real positive numbers with |ℒ|=L<∞. This
set will parameterize the bandwidths  and  for X_C|x
and Y_C|y, respectively, resulting in L estimators in
the ensemble. In other words, we set (l)=l_x^-β and
(l)=l_y^-α with l∈ℒ. While different
parameter sets for  and  can be chosen, we only use one
set here for simplicity of exposition. To achieve the parametric rate,
we need to ensure that the final terms in (<ref>)
are O(1/√(N)). Thus we require the following conditions to
be met: 

    sα   ≥1/2,
    
    sβ   ≥1/2,
    
    1-β d_X-α d_Y   ≥1/2.

For all of these conditions to hold, it is necessary that s≥ d_X+d_Y.
Thus for each estimator in the ensemble we choose (l)=l_x^-1/(2(d_X+d_Y))
and (l)=l_y^-1/(2(d_X+d_Y))
where l∈ℒ. Define w to be a weight vector parameterized
by l∈ℒ with ∑_l∈ℒw(l)=1 and define

    w=∑_l∈ℒw(l)∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/Nh_X_C|x(l),h_Y_C|y(l).

This is the weighted ensemble estimator. From Theorem <ref>,
the bias of w is 

    [ w]     =∑_l∈ℒ∑_i=1^rθ(w(l)l^iN^-i/2(d_X+d_Y))
        +O(√(L)||w||_2(N^-s/2(d_X+d_Y)+N^-1/2)),

where we use θ notation to omit the constants. 

We use the general theory of optimally weighted ensemble estimation
in <cit.> to improve the MSE convergence
rate of the plug-in estimator by choosing the appropriate weights
to cancel the lower order terms in (<ref>): 
 Let ℒ be a set of
real positive numbers with |ℒ|=L<∞ and let J={1,2,…,d_X+d_Y}.
Assume the same conditions in Theorems <ref> and
<ref> hold with (l)=l_x^-1/(2(d_X+d_Y))
and (l)=l_y^-1/(2(d_X+d_Y)).
Assume that s≥ d_X+d_Y and define w as in (<ref>).
Then the MSE of w_0 attains the parametric rate of convergence
of O(1/N) where w_0 is the solution to the following
offline convex optimization problem:

    [                 min_w               ||w||_2;            subject to          ∑_l∈ℒw(l)=1,;                       ∑_l∈ℒw(l)l^i=0, i∈ J. ]



To summarize, if the weights are chosen using eq. (<ref>),
then the weighted ensemble estimator w_0 achieves the parametric
MSE rate. In practice, the optimization problem in (<ref>)
typically results in a very large increase in variance for finite
samples. Thus we use a relaxed version of (<ref>): 

    [                                      min_w                                          ϵ;                                 subject to                               ∑_l∈ℒw(l)=1,;                                            |∑_l∈ℒw(l)l^iN^1/2-i/2(d_X+d_Y)|≤ϵ,  i∈ J,;                                                                         ‖ w‖ _2^2≤ηϵ. ]

The parameter η is chosen to achieve a trade-off between bias
and variance. As shown in <cit.>,
the ensemble estimator w_0 using the resulting weight vector
from the optimization problem in (<ref>) still achieves
the parametric MSE convergence rate under the same assumptions as
described previously. We denote this estimator as GENIE. Algorithm <ref>
summarizes the estimator GENIE.


A similar approach can be used to derive an ensemble estimator for
the case when  and  are purely continuous. Furthermore,
we can define ensemble estimators for both the continuous and the
mixed cases that achieve the parametric MSE rate if s>(d_X+d_Y)/2,
although the optimization problem is more complicated. See Appendix <ref>
for details.

The weights obtained in (<ref>) are optimal in two senses.
First, they are the optimal solution to the problem in (<ref>).
This contrasts with other popular ensemble methods such as random
forests, where the ensemble of learners are equally weighted, and
AdaBoost, where the weights are assigned to different regions of the
feature space based on the training data. The weights are also optimal
in an asymptotic sense. It can be shown that the variance of the ensemble
estimator is bounded by a multiple of ||w||_2 <cit.>.
By minimizing the norm of the weights (or an upper bound on it), we
choose a weight vector that reduces the bias (due to the constraints)
while controlling the variance. Thus the weights are also optimal
in the sense that the bias is reduced to the parametric rate while
the variance is controlled as much as possible given the information
that we have. Since the parametric rate is minimax optimal, this is
also asymptotically optimal for sufficiently smooth densities.

We note that the ensemble estimation approach given here can be compared
to the Jackknife bias correction method <cit.>.
Both approaches use a linear combination of estimators to obtain a
less-biased estimator. However, the standard Jackknife approach uses uniform
weights for the linear combination while the ensemble approach presented
here obtains weights from an optimization problem. This results in
a more computationally efficient procedure as only L estimators
are required for the ensemble approach where L is on the order
of 30-50. The weights can also be computed offline and so solving
the optimization problem contributes little to the total computation
time. In contrast, the standard Jackknife approach requires N different
estimators which is less computationally efficient. 

A more general Jackknife approach such as that in <cit.> shares more similarities with our ensemble method. In this particular work, the authors similarly compute the weights based on an asymptotic bias expansion. However, they do not control the variance via the norm of the weights as we do. Additionally, the Jackknife approach uses a linear combination of estimators with different samples sizes while we use estimators with different bandwidths. Finally, this general Jackknife approach is still more computationally intensive than our ensemble method which computes the weights offline.

At first glance, the weighted ensemble approach discussed in this
section appears to be quite similar to the optimal kernel approaches
used in <cit.>.
However, the weighted ensemble estimation theory we use is applied
to an ensemble of MI estimators after plugging in an ensemble of KDEs
with different bandwidths. So in some sense, we are optimizing the
ensemble of kernels (whose shape is determined by the bandwidth and
the fixed kernel) for the MI estimation problem. In contrast, the
optimal KDE approach first optimizes the kernel for the KDE problem,
and then plugs in the optimized KDE for MI estimation. It is possible
that a proper modification of the ensemble estimation theory could
be applied to a KDE to obtain an optimal KDE and unify these approaches.
This extension is left for future work.



 §.§ Parameter Selection


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


  * Select the minimum and maximum bandwidth parameter to produce density
estimates that satisfy the following: first the minimum bandwidth
should not lead to a zero-valued density estimate at any sample point;
second the maximum bandwidth should be smaller than the diameter of
the support. 

  * Ensure the bandwidths are sufficiently distinct. Similar bandwidth
values lead to a negligible decrease in the bias and many bandwidth
values may increase ||w_0||_2 resulting in an increase in variance
<cit.>.

  * Select L=|ℒ|>|J|=I to obtain a feasible solution for
the optimization problems in (<ref>) and (<ref>).
We find that choosing a value of 30≤ L≤60, and setting ℒ
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-of-thumb given above
are followed. Moon et al. <cit.> 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. In particular, the minimum
and maximum bandwidth parameters can be efficiently selected based
on the k nearest neighbor distances of all data points.

Since the optimal weight w_0 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. <cit.>) can be used to reduce the computational
burden. 



 §.§ Central Limit Theorem


We finish this section with central limit theorems
for the ensemble estimators. This enables us to perform hypothesis
testing on the MI measure.
 Let w^cont be a weighted
KDE ensemble estimator of I_ν(;) when  and  are
continuous with bandwidths h_X(l) and h_Y(l) for each estimator
in the ensemble. Assume that the shaping function g is Lipschitz in both
arguments with Lipschitz constant C_g and that h_X(l), h_Y(l)→0,
N→∞, and Nh_X^d_X(l), Nh_Y^d_Y(l)→∞
for each l∈ℒ. Then for fixed ℒ, and if
𝐒 is a standard normal random variable, 

    (( w^cont-[ w^cont])/√([ w^cont])≤ t)→(𝐒≤ t).


The proof is based on an application of Slutsky's Theorem preceded
by an application of the Efron-Stein inequality (see Appendix <ref>).

For the mixed component case, if 𝒮_X and 𝒮_Y
are finite, then the corresponding ensemble estimators also obey a
central limit theorem. The proof follows by an application of Slutsky's
Theorem combined with Theorem <ref>.
 Let w be a weighted KDE ensemble
estimator of I(;) when  and  contain both continuous
and discrete components. Let the bandwidths for the conditional estimators
be (l) and (l) for each estimator in the ensemble. Assume
that the shaping function g is Lipschitz in both arguments and that , →0,
N→∞, and Nh_X^d_X, Nh_X|y^d_X→∞
for each l∈ℒ and ∀(x,y)∈𝒮_X_D×𝒮_Y_D
with |𝒮_X_D|,|𝒮_Y_D|<∞.
Then for fixed ℒ, 

    (( w-[ w])/√([ w])≤ t)→(𝐒≤ t).





§ APPLICATIONS






 §.§ Simulations


In this section, we validate our theory by estimating the Rényi-α
MI integral (i.e. g(x)=x^α in (<ref>); see <cit.>)
where 𝐗 is a mixture of truncated Gaussian random variables
restricted to the unit cube and 𝐘 is a categorical random
variable that indicates the corresponding truncated Gaussian random
variable that  is drawn from in the mixture. In this setting,
 can be viewed as a classification variable and  contains
the chosen features, which are all continuous in this case. Since
 is purely continuous and  is purely discrete, the MI integral
reduces to the following: 

    I(;)=∑_y∈ S_Yf_Y_D(y)∫(f_X_C(x_C)/f_X_C|Y_D(x_C|y))^αf_X_C|Y_D(x_C|y)dx_C.


We illustrate with Rényi MI as it has received recent interest and the
estimation problem does not reduce to entropy estimation, in contrast
to Shannon MI. Thus this is a clear case where there are no other
nonparametric estimators that are known to achieve the parametric
MSE rate. In fact, to the best of our knowledge, there are no other
nonparametric estimators of Rényi MI that are known to be consistent
in this mixed setting.

We consider two cases. In the first case, 𝐘 has three
possible outcomes (i.e. |𝒮_Y|=3) and respective probabilities
(𝐘=0)=(𝐘=1)=2/5 and (𝐘=2)=1/5.
The conditional covariance matrices are all 0.1× I_d and
the conditional means are, respectively, μ̅_0=0.25×1̅_d,
μ̅_1=0.75×1̅_d, and μ̅_2=0.5×1̅_d,
where I_d 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 ℒ 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). 

Figure <ref> shows the MSE (200 trials) of the plug-in
KDE estimator of the MI integral using a uniform kernel and the optimally
weighted ensemble estimator GENIE for various sample sizes
and for d=4, 6, 9, respectively. The ensemble estimator GENIE
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, 𝐘 has six possible outcomes (i.e.
|𝒮_Y|=6) and respective probabilities (=0)=0.35,
(=1)=0.2, (=2)=(=3)=0.15, (=4)=0.1, and
(=5)=0.05. We chose α=0.5 and d=6. The conditional
covariance matrices are again 0.1× I_d and the conditional
means are, respectively, μ̅_0=0.25×1̅_d, μ̅_1=0.75×1̅_d,
and μ̅_2=0.5×1̅_d, μ̅_3=(0.25×1̅_4^T,0.5×1̅_2^T)^T,
μ̅_4=(0.75×1̅_2^T,0.375×1̅_4^T)^T,
and μ̅_5=(0.5×1̅_4^T,0.25×1̅_2^T)^T.
The results are again given in Figure <ref>. The parameters
for the ensemble estimator and the KDE plug-in estimators are the
same as in the other three plots in Figure <ref>. The
ensemble estimator also outperforms the plug-in estimator in this
setting. 

The estimated negative slopes of the log-log plots in Figure <ref>
are given in Table <ref>. In all settings, both the plug-in
and ensemble estimators outperform their theoretical rates in this
finite-sample regime. However, the rates are generally approaching
the theoretical rates as the dimension increases. It is also clear
from these slopes that the ensemble estimators greatly outperform
the plug-in estimators. We expect the rates to converge to the theoretical
rates as the sample size increases.





 §.§ Application to Single-Cell RNA-Sequencing Data


A common application of MI estimation is to measure the strength of
relationships between different variables, especially in a feature
selection setting. Model aggregation, which includes ensemble methods,
for model selection is a classical problem in statistics <cit.>.
Here we use the GENIE estimator on two different single-cell RNA-sequencing
(scRNA-seq) datasets to demonstrate the estimator's utility for feature
selection. 

Information theory has been used previously in many genomics applications <cit.>.
Single-cell RNA-sequencing data is obtained by measuring the RNA expression
levels in individual (i.e. single) cells <cit.>.
Thousands of genes are typically measured in thousands of cells. This
allows the data to capture the heterogeneity of cell types within
a sample, in contrast with bulk RNA-sequencing methods which effectively
measure the average RNA expression levels within a sample. To correct
for undersampling that is present in scRNA-seq data, we first performed
imputation on both datasets <cit.>. 

For these datasets, we estimated two MI measures: the Rényi MI
and DREMI <cit.>. We define the Rényi
MI to be equal to the Rényi divergence between the joint distribution
of  and  and the product of the marginal distributions.
The DREMI score is a weighted MI developed specifically for analyzing
single-cell data <cit.>. See Appendix <ref>
for further details. Note that no other estimator has been defined
for I_DREMI when the dimension of the continuous component or
components are greater than 1.



  §.§.§ Mouse bone marrow data


We applied GENIE to scRNA-seq data measured from developing mouse
bone marrow cells <cit.>. Estimating mutual
information is commonly done in feature selection where features (in
this case the expression levels of genes) are selected based on the
estimated mutual information between the features  (in this case
the gene expression levels) and the response variable  (in this
case the cell type classification). Features with higher MI are chosen
as they provide more information about the response variable. After
preprocessing, the data contained 10,738 genes measured in 2,730 cells.
In <cit.>, the authors assigned each of the
cells to one of 19 different cell types based on its gene expression
profile. Examples of cell types in this data include erythrocytes,
basophils, and monocytes. 

For this data, we estimated the two different MI measures between
the cell type classification (discrete) and selected groups of genes
(continuous). We estimated the MI for different combinations of genes
selected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways
associated with the hematopoietic cell lineage <cit.>.
Each of these collections contained 8-10 genes. Since the number of
cell types is discrete and the gene expression levels are continuous,
the estimation problem corresponds to estimating the MI between 
and  for the case where  is discrete and  is continuous.
In this problem, |𝒮_Y|=19 and d_X is the number
of genes in the chosen collection.

Table <ref> gives the results. The mean and standard deviation
of the estimated MI (calculated from 1000 bootstrap samples) are reported
for each gene collection including all genes from the four selected
KEGG pathways. Note that the scores for DREMI and Rényi MI are
not directly comparable due to different scaling. The estimated Rényi
MI for these collections is higher than when selecting 8 genes at
random. This is corroborated by classification accuracies obtained
using either a linear SVM classifier or random forests: the classification
accuracies using the KEGG pathways genes are significantly higher
than those obtained using a random set of genes. This suggests the
genes in KEGG pathways associated with the hematopoietic lineage do
provide some information about cell type in this data. Additionally,
the combined genes from all four pathways have the largest estimated
MI for both measures and classification accuracy, which is expected
as genes from different pathways contain information about different
cell types and are thus necessary for distinguishing between cell
types. 

In general, the estimated DREMI when using the KEGG pathways is higher
than the estimated DREMI obtained using random genes. However, several
of these scores are within a standard deviation of the score obtained
from the random genes. Of the four KEGG pathways collections, the
Erythrocyte pathway genes has the largest estimated Rényi MI
and smallest estimated DREMI. Yet, the classification accuracy is
essentially the same as that of the Platelets pathway geneset. These
results highlight the different use cases of these two MI measures.
The Erythrocyte cells are the largest group, containing 1,095 cells.
This suggests that the estimated Rényi MI is biased high for
features relevant for overrepresented groups. In contrast, the DREMI
score appears to be biased low in this case. These results indicate
that the DREMI score may be more appropriate than the Rényi MI
when analyzing less common populations. On the other hand, when less
common populations are not relevant to the analysis, DREMI may not
be as appropriate as other MI measures. These different use cases
highlight the utility of the GENIE estimator in estimating different
MI measures.





  §.§.§ Human embryoid body data


We applied GENIE to scRNA-seq data measured from human embryoid bodies
(EB) collected over a 27-day time course <cit.>. Cells
were sampled at 3-day intervals and then pooled resulting in 5 different
sample collections over time. Thus sample 1 contains cells from days
0 and 3, sample 2 contains cells from days 6 and 9, etc. After preprocessing,
the data contained 17,580 genes measured in 16,825 cells, with each
of the five time samples containing about 2,400 to 4,100 cells. In <cit.>,
the authors identified and analyzed several branches of cells. We
used GENIE to identify genes associated with a neural progenitor (NP)
branch and a neural crest (NC) branch by estimating the Rényi
MI and the DREMI score between the gene expression levels of the cells
in each branch () and the timecourse variable (). This again
corresponds to the case where  is discrete and  is continuous.
For this problem, |𝒮_Y|=5 and d_X is allowed to
vary as described below. Figure <ref> shows PHATE visualizations
of the data highlighted by time sample, and the two branches.


We performed three experiments with each of the branches. For all
experiments, we limited ourselves to genes that are on average nondecreasing
in the branch as time goes on. Thus in each branch, we only considered
the genes such that the correlation between the gene expression level
and time is greater than zero. 

For the first experiment, we estimated the MI scores between the time
course variable and a single gene for all genes in the data (i.e.,
d_X=1). Table <ref> contains the estimated MI scores
of the top 10 genes for each of the measures and branches. Several
of these genes are known to be associated with their respective tissues.
For example, CX3CL1 is often expressed in the brain <cit.>,
SEPT6 has been found to be important for the developing neural tube
in zebra fish <cit.>, SREBF2 is necessary for normal
brain development in mice <cit.>, NR2E1 is predominantly
expressed in the developing brain <cit.>, and ZNF804A
may help regulate early brain development <cit.>.
For the NC branch, multiple HOX genes are listed as having high Rény
MI, all of which are known to be important in the NC <cit.>.
Additionally, RBP1 has been found in enteric nerve NC cells <cit.>,
SHC4 is involved in melanocyte (an NC derivative) development <cit.>,
and PRAME is involved in further differentiation of NC cells <cit.>. 


For comparison, we also used the sure independence screening (SIS)
approach described in <cit.>. This approach reduces
to selecting the genes with the largest correlation with . Table <ref>
shows the top 10 genes for each of the branches and the corresponding
correlation coefficient. Note that only 1/10 of the SIS-selected genes
match with the Rényi MI-selected genes in the NP branch and only
3/10 in the NC branch. None of the DREMI-selected genes match the
SIS-selected genes. Since the SIS approach focuses on linear relationships,
this suggests that our MI estimator is able to effectively detect
strong relationships that are not strictly linear.

None of the DREMI-selected genes match the Rényi MI-selected
genes in both branches. Visualizing the gene expression levels of
the selected genes using PHATE indicates that genes with high DREMI
scores tend to be more localized to a branch while genes with high
Rényi MI may be spread out more (see Figure <ref>
for some examples). This suggests that the DREMI score may be better
than the Rényi MI when the goal is to identify genes that are
uniquely expressed in specific branches. Again, these different use
cases highlight the utility of the GENIE estimator in estimating different
MI measures.


For the second experiment, we used a greedy forward-selection approach
with the GENIE estimator to identify relevant genes for the two branches.
For the third experiment, we used the same forward-selection approach
as in the second experiment except we started by including three or
four relevant genes identified in <cit.> per branch.
In these experiments, we were able to identify several genes with
known associations with these cell types. See Appendix <ref>
for details.

Our results here indicate that GENIE can be useful in identifying
relevant features under multiple settings, even when the variables
are not purely continuous or purely discrete. In particular, since
GENIE accurately identifies previously known gene relationships, we
propose that GENIE can be used to identify unknown gene relationships
for biological discovery. This use can also be extended to other domains
for scientific discovery. 



§ CONCLUSION


We derived the MSE convergence rates for general plug-in KDE-based
estimators of general MI measures between 𝐗 and 𝐘
when they have only continuous components and for the case where 
and/or  contain a mixture of discrete and continuous components.
Using these rates, we defined an ensemble estimator GENIE that achieves
an MSE rate of O(1/N) when the densities are sufficiently smooth.
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.  and  are not both purely
discrete or purely continuous). We also derived a minimax lower bound
on the convergence rate for estimating MI in the continuous case,
derived the asymptotic distribution of the estimator, validated the
superior convergence rate of the ensemble estimator via experiments,
and applied the estimator to analyze feature relevance in single cell
data. We show that the ensemble estimators for the continuous case
achieve the minimax rate for sufficiently smooth densities. Future
work includes extending this approach to k-nn based estimators
which are generally computationally easier than KDE estimators and
extending the minimax rate for the mixed case considered here. We
conjecture that the minimax rates in the mixed case are at least as
slow as those for the continuous case as the mixed case can be decomposed
into a random sum of continuous MI estimators.





ieeetr









§ EXTENDED GENOMICS DETAILS AND RESULTS


Here we provide further details on the genomics experiments. 



 §.§ Imputation and DREMI


Single-cell RNA-sequencing data typically suffers
from undersampling. Therefore, we perform imputation on both datasets
using MAGIC <cit.> prior to estimating the MI
measures. 

Typically, MI measures are weighted by the joint probability density
of  and . In DREMI, the measure is instead weighted by the
conditional probability density of |. This allows DREMI to
measure the strength of the relationship between and  regardless
of differences in population density that often arise in single-cell
data. Since  is continuous and  is discrete for both genomics
applications, DREMI can be defined mathematically as


    I_DREMI(;)     =∑_y∈ S_Y∫ f_Y_D|X_C(y|x)log(f_X_DY_C(x_C,y)/f_X_C(x_C)f_Y_D(y))dx_C
        =∑_y∈ S_Yf_Y_D(y)∫log(f_X_C|Y_D(x_C|y)/f_X_C(x_C))f_X_C|Y_D(x_C|y)/f_X_C(x_C)dx_C.

This measure differs from standard Shannon MI with the inclusion of
the weight 1/f_X_C(x_C) within the integral. While
this does not fit our standard definition of a generalized MI, our
estimation approach allows us to include the inverse of the KDE of
f_X_C when estimating the integral. The proof techniques are
unaffected and therefore our theoretical results still hold.



 §.§ EB Data Extended Results


For the second experiment, we used a greedy forward-selection
approach with the GENIE estimator to identify relevant genes. We first
selected the gene with the highest estimated MI in a given branch
(d_X=1). We then identified the gene that gave the largest MI
when included with the first gene (d_X=2). We then repeated this
to obtain the top 10 genes. The results are shown in Table <ref>.
Rényi MI should never decrease as we add more genes, and we indeed
see this in Table <ref>. Thus the relative increase in
estimated Rényi MI can be used as a measure of the amount of
information each gene adds. Note that for both branches, the largest
increase in Rényi MI occurs within the first four genes and the
inclusion of each subsequent gene adds a decreasing amount of Rényi
MI. However, several of these genes have known associations with their
respective branches. Mutations of HFE are associated with neurological
disorders <cit.> while DOC2A is mainly expressed in
the brain <cit.>. For the NC branch, RGR
is associated with eye development which comes partially from the
neural crest, ITPKB is associated with neurulation <cit.>,
and DPYSL4 is associated with the development of the nervous system <cit.>. 


While the Rényi MI does not decrease with the addition of genes,
DREMI may decrease due to the reweighting caused by using the conditional
distribution instead of the joint. Thus the change in score when adding
genes is less informative for DREMI. For a fixed dimension, however,
the relative DREMI scores are informative and thus can be used to
identify relevant genes using the forward-selection approach. Using
this approach with DREMI, we identified several genes with known associations
such as HFE (also identified with Rényi MI) and BRD9 <cit.>
with the NP branch, and CFL1 <cit.> and BMP8B <cit.>
with the NC branch.

We also performed a forward-selection variant on SIS. We first selected
the gene with the highest SIS score (correlation coefficient in this
case). We then performed regression with this gene and the time course
variable . We then calculated the SIS score between all of the
other genes individually and the regression residuals to select the
next gene. This process was repeated to obtain a list of the top ten
genes in Table <ref>. Since the SIS criteria is scale-invariant,
this can sometimes result in an increase in the correlation coefficient
as more genes are included, although generally we expect the correlation
to decrease. Thus it is somewhat difficult to assess using SIS the
amount of information added by including each gene. In this case,
the MI and SIS approaches identified unique genes with no shared overlap
in either branch, again suggesting that our MI approaches are identifying
nonlinear relationships. 

For the third experiment, we used the same forward-selection approach
as in the second experiment except we started by including three or
four relevant genes identified in <cit.>. These genes
were NKX2-8, EN2, and SOX1 for the NP branch, and PAX3, FOXD3, SOX9,
and SOX10 for the NC branch. The results are presented in Table <ref>.
Interestingly, including these “preset” genes results in a larger
overall Rényi MI and DREMI in the NP branch than when using a
purely greedy approach (Table <ref>) while the opposite
is true for the NC branch. Additionally, the identified genes are
all different from the purely greedy approach. However, many of them
are known to be associated with their respective tissues. PTK6 affects
neurite extension <cit.>, SIRT3 regulates
mitochondria in the brain during development <cit.>,
HEY1 is expressed in neural precursor cells <cit.>,
BMPR1B is important for brain development <cit.>,
FAM72C is enriched in cortical neural progenitors <cit.>,
PPFIA4 is involved in neural development <cit.>,
LRSAM1 is related to enteric NC cells <cit.>,
LINC00327 is associated with regulating neuroblasts <cit.>,
and GRHPR is associated with human eye development <cit.>.





§ MI ENSEMBLE ESTIMATION EXTENSIONS


In this appendix, we present several extensions
of the ensemble estimation approach to MI. First, we show how to apply
the theory to the purely continuous case. We then show how the theory
can be applied to obtain estimators that achieve the parametric rate
for less smooth densities.



 §.§ Continuous Random Varables


We can apply Theorem 3 in <cit.>
to obtain a version of the GENIE MI estimator that achieves the parametric
rate for the case when 𝐗 and 𝐘 are purely
continuous. For completeness, we repeat the theorem and its proof here. For a general
estimation problem, let N be the number of available samples and
let ℒ={ l_1,…,l_L} be a set of
index values. For an indexed ensemble of estimators {𝐄̂_l} _l∈ℒ
of a parameter E, the weighted ensemble estimator with weights
w={ w(l_1),…,w(l_L)}
satisfying ∑_l∈ℒw(l)=1 is defined as 

    𝐄̂_w=∑_l∈ℒw(l)𝐄̂_l.

 Consider the following conditions on {𝐄̂_l} _l∈ℒ: 


  * 𝒞.1 The bias is expressible as 

    [𝐄̂_l]=∑_i∈ Jc_iψ_i(l)ϕ_i,d(N)+O(1/√(N)),

where c_i are constants depending on the underlying density and
are independent of N and l, J={ i_1,…,i_I}
is a finite index set with I<L, and ψ_i(l) are basis functions
depending only on the parameter l and not on the sample size N. 

  * 𝒞.2 The variance is expressible as 

    [𝐄̂_l]=c_v(1/N)+o(1/N).


Assume conditions 𝒞.1
and 𝒞.2 hold for an ensemble of estimators {𝐄̂_l} _l∈ℒ.
Then there exists a weight vector w_0 such that the MSE of the
weighted ensemble estimator attains the parametric rate of convergence:

    𝔼[(𝐄̂_w_0-E)^2]=O(1/N).

The weight vector w_0 is the solution to the following convex
optimization problem:

    [                           min_w                         ||w||_2;                      subject to                    ∑_l∈ℒw(l)=1,;                                 γ_w(i)=∑_l∈ℒw(l)ψ_i(l)=0, i∈ J. ]




Due to condition 𝒞.1, the bias of the weighted ensemble estimator is

    [𝐄̂_w]=∑_i∈ Jc_iγ_w(i)ϕ_i,d(N)+O(√(L)||w||_2/√(N). )

Denote the covariance matrix of the ensemble of estimators as Σ_L. By the Cauchy-Schwarz inequality and condition 𝒞.2, the entries of Σ_L are O(1/N). The variance of the weighted estimator is then bounded above as 

    [𝐄̂_w]=w^TΣ_L w≤Trace(Σ_L N)||w||_2^2/N=c_v L||w||_2^2/N+o(1/N).


The optimization problem in (<ref>) zeroes out the lower-order bias terms and limits the ℓ_2 norm of the weight vector w to prevent the variance from exploding. This results in an MSE of O(1/N) as long as the dimension d is fixed and L is fixed and independent of the sample size N. A solution is guaranteed to (<ref>) as long as L>I and the vectors a_i=[ψ_i(l_1),…,ψ_i(l_L)] are linearly independent.



As before, (<ref>) typically results in an ensemble
estimator with a large variance. We can relax this optimization problem
and obtain an estimator that still obtains the parametric rate: 

    [                           min_w                               ϵ;                      subject to                    ∑_l∈ℒw(l)=1,;                                 |γ_w(i)N^1/2ϕ_i,d(N)|≤ϵ,  i∈ J,;                                                   ‖ w‖ _2^2≤ηϵ. ]


We can use (<ref>) to obtain a GENIE estimator for the
purely continuous case. Theorem <ref> indicates that we
need h_X^d_Xh_Y^d_Y∝ N^-1/2 for the O(1/(Nh_X^d_Xh_Y^d_Y))
terms to be O(1/√(N)). We consider the more general case where
the parameters may differ for h_X and h_Y. Let ℒ_X
and ℒ_Y be sets of real, positive numbers with |ℒ_X|=L_X
and |ℒ_Y|=L_Y. For each estimator in the ensemble,
choose l_X∈ℒ_𝒳 and l_Y∈ℒ_Y and
set h_X(l_X)=l_XN^-1/(2(d_X+d_Y)) and h_Y(l_Y)=l_YN^-1/(2(d_X+d_Y)).
Define the matrix w s.t. ∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)=1.
From Theorems <ref> and <ref>, conditions
𝒞.1 and 𝒞.2 are satisfied if s≥ d_X+d_Y
with ψ_i,j(l_X,l_Y)=l_X^il_Y^j and ϕ_i,j(N)=N^-(i+j)/(2(d_X+d_Y))
for 0≤ i,j≤ d_X+d_Y s.t. 0<i+j≤ d_X+d_Y. The
optimal weight w_0 is calculated using (<ref>).
The resulting estimator 

    w_0^cont=∑_l_X∈ℒ_X,l_Y∈ℒ_Yw_0(l_X,l_Y)h_X(l_X),h_Y(l_Y)

achieves the parametric MSE rate when s≥ d_X+d_Y. We denote
this estimator as GENIE^cont.



 §.§ Less Smooth Densities


The GENIE estimators GENIE and GENIE^cont
are guaranteed to achieve the parametric convergence rate as long
as s≥ d_X+d_Y. Here we derive ensemble estimators of MI
that achieve the parametric rate under less strict smoothness assumptions
on the densities. To derive the ensemble estimators for less smooth
densities, we need a different expansion of the bias that includes
some higher order terms. We present those results here for the continuous
and mixed cases and show how to apply the ensemble estimation theory
to obtain the parametric MSE convergence rate when the densities are
less smooth (s>(d_X+d_Y)/2).



  §.§.§ Continuous Random Variables


We first consider the case where  and  are both purely continuous.
Consider the following result on the bias of the plug-in estimator: 
Assume that assumptions
𝒜.0-𝒜.5 hold. Then for λ≥2 a
positive integer, the bias of  is 

    []     =    ∑_m,n=0
    
    i+j+m+n≠0
    ^λ∑_i,j=0^rc_11,i,j,m,nh_X^ih_Y^j/(Nh_X^d_X)^m(Nh_Y^d_Y)^n
           +∑_m=1^λ∑_i=0^r∑_j=0^rc_13,i,j,mh_X^ih_Y^j/(Nh_X^d_Xh_Y^d_Y)^m
           +O(h_X^s+h_Y^s+1/(Nh_X^d_Xh_Y^d_Y)^λ).


The proof is given in Appendix <ref>. Note that no
extra assumptions are required to achieve this result compared to
Theorem <ref>. However, we elected to retain these results
for the appendix to simplify the presentation in the main paper.

We now use these results to define a new ensemble estimator. Set δ>0
and let ℒ_X and ℒ_Y be sets of real,
positive numbers with |ℒ_X|=L_X and |ℒ_Y|=L_Y.
For each estimator in the ensemble, choose l_X∈ℒ_𝒳
and l_Y∈ℒ_Y and set h_X(l_X)=l_XN^-1/(d_X+d_Y+δ)
and h_Y(l_Y)=l_YN^-1/(d_X+d_Y+δ). Then conditions
𝒞.1 and 𝒞.2 are satisfied if s≥(d_X+d_Y+δ)/2
and λ≥(d_X+d_Y+δ)/δ with ψ_1,i,j,m,n(l_X,l_Y)=l_X^i-md_Xl_Y^j-nd_Y
and ϕ_1,i,j,m,n(N)=N^-i+j+m(d_Y+δ)+n(d_X+δ)/d_X+d_Y+δ
for 0<i+j+m(d_Y+δ)+n(d_X+δ)≤d_X+d_Y+δ/2
and the terms ψ_2,i,j,m(l_X,l_Y)=l_X^i-md_Xl_Y^j-md_Y
and ϕ_2,i,j,m(N)=N^-i+j+mδ/d_X+d_Y+δ
for m≥1 and i+j+mδ≤d_X+d_Y+δ/2.
These all correspond to terms that converge to zero slower than N^-1/2
when left uncorrected. The optimal weight w_0 is again calculated
using (<ref>) and the resulting ensemble estimator achieves
the parametric MSE convergence rate when s≥(d_X+d_Y+δ)/2.
Since δ can be chosen arbitrarily close to zero, the parametric
rate can be achieved theoretically as long as s>(d_X+d_Y)/2. 



  §.§.§ Mixed Random Variables


We now consider the case where 𝐗 and 𝐘 may
have any mixture of continuous and discrete components. We have a
similar result on the bias as in Theorem <ref>.
Here we assume that =l_X𝐍_x^-β and =l_Y_y^-α
with 0<β<1/d_X, 0<α<1/d_Y, and
l_X,l_Y>0.
Assume that the same assumptions hold as in
Theorem <ref>. Then for λ≥2 a positive
integer, the bias of h_X_C|X_D,h_Y_C|Y_D is

    [h_X_C|X_D,h_Y_C|Y_D]     =∑_m,n=0
    
    i+j+m+n≠0
    ^λ∑_i,j=0^rc_14,i,j,m,nl_X^il_Y^jN^-iβ-jα/(l_X^d_XN^1-β d_X)^m(l_Y^d_YN^1-α d_Y)^n
        +∑_m=1^λ∑_i=0^r∑_j=0^rc_15,i,j,ml_X^il_Y^jN^-iβ-jα/(l_X^d_Xl_Y^d_YN^1-β d_X-α d_Y)^m
        +O(N^-sβ+N^-sα+1/(N^1-β d_X-α d_Y)^λ).


The proof is given in Appendix <ref>.

We now use these results to define a new ensemble estimator in the
mixed case. The procedure is similar to the continuous case. Set δ>0
and let ℒ_X and ℒ_Y be sets of real,
positive numbers with |ℒ_X|=L_X and |ℒ_Y|=L_Y.
For each estimator in the ensemble, choose l_X∈ℒ_𝒳
and l_Y∈ℒ_Y and set (l_X)=l_X𝐍_x^-1/(d_X+d_Y+δ)
and (l_Y)=l_Y_y^-1/(d_X+d_Y+δ). Conditions
𝒞.1 and 𝒞.2 are satisfied if s≥(d_X+d_Y+δ)/2
and λ≥(d_X+d_Y+δ)/δ. The first set of terms
in the optimization problem are ψ_1,i,j,m,n(l_X,l_Y)=l_X^i-md_Xl_Y^j-nd_Y
and ϕ_1,i,j,m,n(N)=N^-i+j+m(d_Y+δ)+n(d_X+δ)/d_X+d_Y+δ
for 0<i+j+m(d_Y+δ)+n(d_X+δ)≤d_X+d_Y+δ/2.
The second set of terms are ψ_2,i,j,m(l_X,l_Y)=l_X^i-md_Xl_Y^j-md_Y
and ϕ_2,i,j,m(N)=N^-i+j+mδ/d_X+d_Y+δ
for m≥1 and i+j+mδ≤d_X+d_Y+δ/2.
The optimal weight w_0 is again calculated using (<ref>)
and the resulting ensemble estimator achieves the parametric MSE convergence
rate when s≥(d_X+d_Y+δ)/2. Since δ can be chosen
arbitrarily close to zero, the parametric rate can be achieved theoretically
as long as s>(d_X+d_Y)/2. 

The modified estimators defined in this section have better statistical
properties than the original GENIE estimators defined in Section <ref>
and Appendix <ref> 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 (<ref>) is larger for the modified estimator.
In theory, this could lead to larger variance although this is not
necessarily true in practice according to divergence estimation experiments
in <cit.>.



§ THE BOUNDARY CONDITION


Here we prove that under certain smoothness assumptions on a kernel
with either rectangular or circular support, the boundary assumption 𝒜.5 is
satisfied for densities with the unit cube as its support set. We will
prove the following, more general result:

 Let K be a kernel function  with rectangular or circular support,
i.e., K(u)=0 for ‖ u‖_1 >1 or ‖ u‖_2 >1, respectively. Let K∈Σ(s,H_K) in the interior of its support. Let p_x(u):ℝ^d→ℝ
be a polynomial in u of order |q|≤ r=⌊ s⌋
whose coefficients are a function of x and are r-|q| times differentiable.
Then for 𝒮=[0,1]^d and any positive integers t and
m, we have that

    ∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮K^l(u)p_x(u)du)^tdx=v_t(h),

 where v_t(h) admits the expansion 

    v_t(h)=∑_i=1^r-|q|e_i,q,t,lh^i+o(h^r-q),

 for some constants e_i,q,t,l.


Before proving this, we will relate this result to assumption 𝒜.5.


Assumption 𝒜.5 is satisfied under the same conditions as in Theorem <ref>, i.e., K∈Σ(s,H_K) is a kernel function with either rectangular or circular support and 𝒮=[0,1]^d.

This follows immediately from assumption 𝒜.2 that f is in the Hölder class Σ(s,H), which implies that D^q f(x) is r-|q| times differentiable. Thus equation (<ref>) implies 𝒜.5.





 §.§ Proof of Theorem <ref>: Rectangular Support Kernels

We will first consider points that are boundary points due to a single coordinate and then extend to the general case where multiple coordinates are close to the boundary.



  
*Single Coordinate Boundary Point Since K∈Σ(s,H_K), we can take a Taylor series expansion of K^l(u) around zero to obtain

    K^l(u)=p_K,l(u)+o(||u||_2^r),

where p_K,l(u) is a  polynomial function of u with degree r.


Consider points x that are boundary points by virtue
of a single coordinate x_i such that x_i+u_ih∉.
Without loss of generality, assume that x_i+u_ih>1. After performing the above Taylor series expansion of K^l, the inner
integral in (<ref>) can then be evaluated first with respect to
all coordinates other than i. Since all of these coordinates lie
within the support, the inner integral over these coordinates will
amount to integration of the polynomial p_x(u) over a symmetric
d-1 dimensional rectangular region |u_j|≤ 1
for all j≠ i. This yields a function ∑_m=1^q+rp̃_m(x)u_i^m+o(u_i^r+1)
where the coefficients p̃_m(x) are each r-|q| times differentiable
wrt x.

With respect to the u_i coordinate, the inner integral will have
limits from 1-x_i/h to 1 for some 1>x_i>1-h.
Consider the p̃_m(x)u_i^m monomial term. The inner
integral wrt this term yields (ignoring the o(·) term for now)

    ∑_m=1^|q|+rp̃_m(x)∫_1-x_i/h^1u_i^mdu_i=∑_m=1^|q|+rp̃_m(x)1/m+1(1-(1-x_i/h)^m+1).

Raising the right hand side of (<ref>) to the power of t
results in an expression of the form 

    ∑_j=0^(|q|+r)tp̌_j(x)(1-x_i/h)^j,

where the coefficients p̌_j(x) are r-|q| times differentiable
wrt x. Integrating (<ref>) over all the coordinates in
x other than x_i results in an expression of the form 

    ∑_j=0^(|q|+r)tp̅_j(x_i)(1-x_i/h)^j,

 where again the coefficients p̅_j(x_i) are r-|q| times
differentiable wrt x_i. Note that since the other cooordinates
of x other than x_i are far away from the boundary, the coefficients
p̅_j(x_i) are independent of h. To evaluate the integral
of (<ref>), consider the r-|q| term Taylor series expansion
of p̅_j(x_i) around x_i=1, where we use a smooth extension of the function and its derivatives to the boundary. This will yield terms
of the form 

    ∫_1-h^1(1-x_i)^j+k/h^kdx_i    =    .-(1-x_i)^j+k+1/h^k(j+k+1)|_x_i=1-h^x_i=1
        =    h^j+1/(j+k+1),

for 0≤ j≤ r-|q|, and 0≤ k≤ (|q|+r)t. By a similar analysis, the o(·) terms from the Taylor expansion of the kernel result in o(h^r). Combining terms results
in the expansion v_t(h)=∑_i=1^r-|q|e_i,q,th^i+o(h^r-|q|).



  
*Multiple Coordinate Boundary Point The case where multiple coordinates of the point x are near the
boundary is a straightforward extension of the single boundary point
case so we only sketch the main ideas here. As an example, consider
the case where 2 of the coordinates are near the boundary. Assume
for notational ease that they are x_1 and x_2 and that x_1+u_1h>1
and x_2+u_2h>1. The inner integral in (<ref>)
can again be evaluated first wrt all coordinates other than 1 and
2. This yields a function ∑_m,j=1^q+rp̃_m,j(x)u_1^mu_2^j
where the coefficients p̃_m,j(x) are each r-q times
differentiable wrt x. Integrating this wrt x_1 and x_2
and then raising the result to the power of t yields a double sum
similar to (<ref>). Integrating this over all the coordinates
in x other than x_1 and x_2 gives a double sum similar
to (<ref>). Then a Taylor series expansion of the coefficients
and integration over x_1 and x_2 yields the result.




 §.§ Proof of Theorem <ref>: Circular Support Kernels

The case where the kernel K has a circular support is more complex than when the support is rectangular. We will again first consider points that are boundary points due to a single coordinate. We will then extend to the general case where multiple coordinates are close to the boundary.



  
*Single Coordinate Boundary Point Consider points x that are
boundary points due to a single coordinate x_i s.t. x_i+u_ih∉𝒮.
Without loss of generality, assume that x_i+u_ih>1. We focus
first on the inner integral in (<ref>). We will use the
following lemma: 

Let D_d(ρ) be a d-sphere with radius
ρ and let ∑_i=1^dn_i=q. Then 

    ∫_D_d(ρ)u_1^n_1u_2^n_2… u_d^n_ddu_1… du_d=Cρ^d+q,

 where C is a constant that depends on the n_is and d.


We convert to d-dimensional spherical coordinates to handle the
integration. Let r be the distance of a point u from the origin.
We nave d-1 angular coordinates ϕ_i where ϕ_d-1
ranges from 0 to 2π and all other ϕ_i range from 0
to π. The conversion from the spherical coordinates to Cartesian
coordinates is then 

    u_1    =     rcos(ϕ_1)
    
    u_2    =     rsin(ϕ_1)cos(ϕ_2)
    
    u_3    =     rsin(ϕ_1)sin(ϕ_2)cos(ϕ_3)
    ⋮
    
    u_d-1    =     rsin(ϕ_1)⋯sin(ϕ_d-2)cos(ϕ_d-1)
    
    u_d    =     rsin(ϕ_1)⋯sin(ϕ_d-2)sin(ϕ_d-1).

 The spherical volume element is then 

    r^d-1sin^d-2(ϕ_1)sin^d-3(ϕ_1)⋯sin(ϕ_d-1)dr dϕ_1 dϕ_2⋯ dϕ_d-1.

 Combining these results gives 

    ∫_D_d(ρ)u_1^n_1u_2^n_2… u_d^n_ddu_1… du_d
        =    ∫_0^ρ∫_o^2π∫_0^π⋯∫_0^πr^q+d-1[sin^q-n_1+d-2(ϕ_1)sin^q-n_1-n_d+d-3(ϕ_2)⋯.
          .sin^n_d+n_d-1+1(ϕ_d-2)sin^n_d(ϕ_d-1)][cos^n_1(ϕ_1)⋯cos^n_d(ϕ_d-1)]dϕ_1⋯ dϕ_d-1dr
        =     Cρ^q+d.



Since K∈Σ(s,H_K), we can take a Taylor series expansion of K^l(u) around zero to obtain

    K^l(u)=p_K,l(u)+o(||u||_2^r),

where p_K,l(u) is a  polynomial function of u with degree r.

The region of integration for the inner integral in (<ref>)
corresponds to a hyperspherical cap with radius 1 and height of
1-x_i/h. The inner integral can be calculated using an
approach similar to that used in <cit.> to calculate
the volume of a hyperspherical cap. It is obtained by integrating
the polynomial p_x(u) over a d-1-sphere with radius sinθ
and height element dcosθ. This is done using Lemma <ref>.
We then integrate over θ which has a range of 0 to ϕ=cos^-1(1-x_i/h). 
Thus we have 

    ∫_u:x+uh∉𝒮K^l(u)p_x(u)du     =    ∫_u:||u||_2≤1,x+uh∉𝒮(p_K,l(u)+o(||u||_2^r))p_x(u)du 
        =    ∑_m=0^|q|+rp̃_m,l(x)∫_0^ϕsin^m+d-1(θ)sinθcos^mθ dθ+o( ∫_u:||u||_2≤1,x+uh∉𝒮||u||_2^rdu)
        =    ∑_m=0^|q|+rp̃_m,l(x)∫_0^ϕsin^m+d(θ)cos^mθ dθ+o( ∫_u:||u||_2≤1,x+uh∉𝒮||u||_2^rdu),

where p̃_m(x) is the polynomial coefficient corresponding to u_i^m after the Taylor series expansion of K^l and integrating over the d-1-sphere. 

We will focus on the first term in Eq. <ref> as the o(·) term will follow similarly as a polynomial function of u. From standard integral tables, we get that for n≥2 and m≥0

    ∫_0^ϕsin^nθcos^mθ dθ=-sin^n-1ϕcos^m+1ϕ/n+m+n-1/n+m∫_0^ϕsin^n-2θcos^mθ dθ.

If n=1, then we get 

    ∫_0^ϕsinθcos^mθ dθ=1/m+1-cos^m+1ϕ/m+1.

Since ϕ=cos^-1(1-x_i/h), we have 

    cosϕ    =    1-x_i/h,
    sinϕ    =    √(1-(1-x_i/h)^2).

Therefore, if n is odd, we obtain 

    ∫_0^ϕsin^nθcos^mθ dθ=∑_ℓ=0^(n-1)/2c_ℓ(√(1-(1-x_i/h)^2))^2ℓ(1-x_i/h)^m+1+c,

where the constants depend on m and n.

If n is even and m>0, then the final term in the recursion in
(<ref>) reduces to 

    ∫_0^ϕcos^mθ dθ=cos^m-1ϕsinϕ/m+m-1/m∫_0^ϕcos^m-2θ dθ.

If m=2, then 

    ∫_0^ϕcos^2θ dθ    =    ϕ/2+1/4sin(2ϕ)
        =    ϕ/2+1/2sinϕcosϕ.

Therefore, if n and m are both even, then this gives 

    ∫_0^ϕsin^nθcos^mθ dθ    =    ∑_ℓ=0^(n-2)/2c_ℓ^'(√(1-(1-x_i/h)^2))^2ℓ+1(1-x_i/h)^m+1+c^'cos^-1(1-x_i/h)
           +∑_ℓ=0^(m-2)/2c_ℓ^”(√(1-(1-x_i/h)^2))(1-x_i/h)^2ℓ+1.

On the other hand, if n is even and m is odd, we get 

    ∫_0^ϕsin^nθcos^mθ dθ    =    ∑_ℓ=0^(n-2)/2c_ℓ^”'(√(1-(1-x_i/h)^2))^2ℓ+1(1-x_i/h)^m+1
           +∑_ℓ=0^(m-1)/2c_ℓ^””(√(1-(1-x_i/h)^2))(1-x_i/h)^2ℓ.


If d is odd, then combining (<ref>) and (<ref>)
with (<ref>) gives 

    ∑_m=0^|q|+rp̃_m,l(x)∫_0^ϕsin^m+d(θ)cos^mθ dθ =  ∑_m=0^|q|+r∑_ℓ=0^d+|q|p_m,ℓ,l(x)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m,

where the coefficients p_m,ℓ,l(x) are r-|q| times differentiable
wrt x. Similarly, if d is even, then 

    ∑_m=0^|q|+rp̃_m,l(x)∫_0^ϕsin^m+d(θ)cos^mθ dθ    =    ∑_m=0^|q|+r∑_ℓ=0^d+|q|p_m,ℓ,l^'(x)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m
           +p^'(x)cos^-1(1-x_i/h),

where again the coefficients p_m,ℓ,l^'(x) and p^'(x)
are r-q times differentiable wrt x. Raising (<ref>)
and (<ref>) to the power of t gives respective expressions
of the form 

    ∑_m=0^(|q|+r)t∑_ℓ=0^(d+|q|)tp̌_m,ℓ,l(x)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m,


    ∑_m=0^(|q|+r)t∑_ℓ=0^(d+|q|)t∑_n=0^tp̌_m,ℓ,l,n(x)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m(cos^-1(1-x_i/h))^n,

where the coefficients p̌_m,ℓ,l(x) and p̌_m,ℓ,l,n(x)
are all r-q times differentiable wrt x. Integrating (<ref>)
and (<ref>) over all the coordinates in x except
for x_i affects only the p̌_m,ℓ,l(x) and p̌_m,ℓ,l,n(x)
coefficients, resulting in respective expressions of the form 

    ∑_m=0^|q|t∑_ℓ=0^(d+|q|)tp̅_m,ℓ,l(x_i)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m,


    ∑_m=0^|q|t∑_ℓ=0^(d+|q|)t∑_n=0^tp̅_m,ℓ,l,n(x_i)(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m(cos^-1(1-x_i/h))^n.


The coefficients p̅_m,ℓ,l(x_i) and p̅_m,ℓ,l,n(x_i)
are r-q times differentiable wrt x_i. Since the other coordinates
of x other than x_i are far away from the boundary, the coefficients
are independent of h. For the integral wrt x_i of (<ref>),
taking a Taylor series expansion of p̅_m,ℓ(x_i) around
x_i=1 (again using a smooth extension of the function and its derivatives to the boundary) yields terms of the form 

    ∫_1-h^1(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m+jh^jdx_i    =     h^j+1∫_0^1(1-y_i)^ℓ/2y_i^m+j-1/2dy_i
        =     h^j+1B(ℓ+2/2,m+j+1/2),

where 0≤ j≤ r-q, 0≤ℓ≤(d+|q|)t, 0≤ m≤ |q|t,
and B(x,y) is the beta function. Note that the first step uses
the substitution of y_i=(1-x_i/h)^2.

If d is even (i.e. (<ref>)), a simple closed-form
expression is not easy to obtain due to the cos^-1(1-x_i/h)
terms. However, by similarly applying a Taylor series expansion to
p̅_m,ℓ,l,n(x_i) and substituting y_i=1-x_i/h
gives terms of the form of 

    ∫_1-h^1(√(1-(1-x_i/h)^2))^ℓ(1-x_i/h)^m+j(cos^-1(1-x_i/h))^nh^jdx_i
        =     h^j+1∫_0^1(1-y_i^2)^ℓ/2y_i^m+j(cos^-1y_i)^ndy_i
        =     h^j+1c_ℓ,m,j,n,

for 0≤ j≤ r-|q|, 0≤ℓ≤(d+|q|)t, 0≤ m≤ (|q|+r)t,
and 0≤ n≤ t. By a similar analysis, the o(·) term in Eq. <ref> results in o(h^r). Combining terms results in the expansion v_t(h)=∑_i=1^r-|q|e_i,q,t,lh^i+o(h^r-|q|).



  
*Multiple Coordinate Boundary Point

The case where multiple coordinates
of the point x are near the boundary is a fairly straightforward
extension of the single boundary point case. Consider the case where
2 of the coordinates are near the boundary, e.g., x_1 and x_2
with x_1+u_1h>1 and x_2+u_2h>1. The region of integration
for the inner integral can be decomposed into two parts: a hyperspherical
cap wrt x_1 and the remaining area (denoted, respectively, as
A_1 and A_2). The remaining area A_2 can be decomposed
further into two other areas: a hyperspherical cap wrt x_2 (denoted
B_1) and a height chosen s.t. B_1 just intersects A_1
on their boundaries. Integrating over the remainder of A_2 is
achieved by integrating along x_2 over d-1-dimensional hyperspherical
caps from the boundary of B_1 to the boundary of A_2. Thus
integrating over these regions yields an expression similar to (<ref>).
Following a similar procedure will then yield the result. 



 §.§ Proof of Proposition <ref>


Recall that we assume that the derivatives of the density f vanish near the boundary of the density support set. We will use the following lemma, which is modified from equation (8) in Singh and Poczos <cit.>:


<cit.> Let the density support set 𝒮 be the unit cube with d≥2 and let the derivatives of the density f vanish near the boundary of the support set, which is denoted as ∂𝒮. Assume that f belongs to the Hölder set of order s, i.e. f∈Σ(H,s). Assume that q with 0≤|q|≤ r and x∈ℬ:={x∈𝒮|dist(x,∂𝒮)≤ h }. Then if h≤ 1/2,

    | D^q f(x)|≤ Hh^s-|q|.



We include a proof here for completeness. We will use an induction argument on |q| as |q| decreases from r to 0. Let y∈∂𝒮 such that ||x-y||≤ h. For |q|=r, we have

    |D^qf(x)|   =|D^qf(x)-D^qf(y)| 
       ≤ H||x-y||^s-|q|
       ≤ Hh^s-|q|,

where the first step comes from the fact that D^qf(y)=0, the second step uses the Hölder condition, and the last step follows from the choice of y. Suppose we have the desired bound for derivatives of order |q|+1. Let x∈ℬ and let u=(0,…,0,± 1,0,…,0)∈ℝ^d, where u_j=± 1 for some j∈[d]. Then we can write x=y+hu for some y∈ℬ. Furthermore, the point y+tu∈ℬ for t∈[0,h]. To see how this is true, note that if x∈ℬ then there exists at least one coordinate x_i such that either 0≤ x_i≤ h or 1-h≤ x_i≤ 1. Pick any other coordinate x_j. Then choose the sign of u_j such that y=x-hu∈𝒮. This is possible since h≤ 1/2.  Since y is constructed by moving from x∈ℬ in parallel with the boundary, then y∈ℬ as well as y+tu∈ℬ.

Based on this construction of y and u, we have

    |D^q f(y+hu)|   ≤∫_0^h | ∂/∂ x_jD^q f(y+tu)|dt 
       ≤∫_0^h Hh^s-(|q|+1)dt 
       =Hh^s-|q|.

The desired result then follows by induction on |q|.




From Lemma <ref>, 

    |∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮K^l(u)u^q D^q f(x) du)^tdx| ≤(Hh^s-|q|||K||_∞)^t| ∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮u^q du)^tdx|.


Consider the case when K has rectangular support from (without loss of generality) -1 to 1 in each dimension. We consider the setting when x is close to the 1 boundary in a single coordinate x_i. The inner integral in (<ref>) then reduces to 

    (C∫_1-x_i/h^1 u_i^m du_i)^t     = (C/m+1(1-(1-x_i/h)^m+1))^t 
        = (C/m+1)^t∑_j=0^t (-1)^jtj(1-x_i/h)^(m+1)j,

where C is a constant that comes from integrating the other terms of u^q and m≤ |q| is the exponent of u_i in u^q. Based on equation (<ref>), taking the integral of this result with respect to x will then yield terms of the form of 
h^j+1 for j≥ 0. The extension to the case with multiple boundary coordinates is similar to that in the proof of Theorem <ref>.

Combining this with the results in equations (<ref>) and (<ref>) gives that 

    |∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮K^l(u)u^q D^q f(x) du)^tdx|=O(h^s-|q|+1)=o(h^r-|q|).

Therefore, assumption 𝒜.5 is satisfied with the coefficients in the expansion being equal to zero.

The proof for the circular support kernel follows similar arguments also adapted from the proof of Theorem <ref>.



 §.§ Assumption 𝒜.5 and the Truncated Gaussian Distribution


Here we show that the truncated Gaussian distribution satisfies assumption 𝒜.5. For simplicity, we will consider the univariate case for s=2 and the uniform kernel. We restrict the standard normal distribution to the interval [-1,1]. We focus on the +1 boundary as the -1 boundary follows a similar procedure. This gives 

    ∫_x∈𝒮(∫_u:K(u)>0, x+uh∉𝒮K^l(u)u^q D^q f(x) du)^tdx    =∫_1-h^1(∫_1-x/h^1 d^q/dx^qf(x)u^q du)^tdx 
        = ∫_1-h^1( (1-(1-x/h)^q+1)d^q/dx^qf(x)/q+1)^tdx.

We can extend the derivatives of the function f to the boundary by simply using the derivatives prior to truncation. By a Taylor series expansion around x=1, we have

    f(x)   =f(1)+f'(1)(1-x)+f”(1)/2(1-x)^2+o((1-x)^2), 
    d/dxf(x)   =f'(1)+f”(1)(1-x)+o(1-x), 
    d^2/dx^2f(x)   =f”(1)+o(1).

Therefore, the integral with respect to x in (<ref>) has terms of the form of (1-x)^j+k/h^k which integrates to h^j+1/j+k+1 (see (<ref>)). Combining these results gives the desired expansion. As an example, suppose t=1 and q=1. Then (<ref>) becomes 

    1/2∫_1-h^1(1-(1-x/h)^2)(f'(1)+f”(1)(1-x)+o(1-x))dx.


Distributing the expression for d/dxf(x) through and evaluating the integrals separately gives

    1/2∫_1-h^1(f'(1)+f”(1)(1-x)+o(1-x))dx=1/2(f'(1)h+f”(1)h^2/2)+o(h^2),


    1/2∫_1-h^1(1-x/h)^2(f'(1)+f”(1)(1-x)+o(1-x))dx=1/2( f'(1)h/3+f”(1)h^2/4)+o(h^2).

Taking the difference between these expressions gives the final expansion:

    f'(1)h/3+f”(1)h^2/8+o(h^2).




§ PROOF OF THEOREM <REF> (CONTINUOUS BIAS)


Here we prove the results shown in Theorem <ref>,
which includes the results for Theorem <ref>. In this section
and throughout all other proofs, let Z denote the conditional
expectation given 𝐙. 

The bias of  is

    []=[g( X(𝐗) Y(𝐘)ν_1ν_2/ Z(𝐗,𝐘)ν_3)-g(f_X(𝐗)f_Y(𝐘)ν_1ν_2/f_XY(𝐗,𝐘)ν_3)],

where 𝐗 and 𝐘 are drawn jointly from f_XY.
We will derive an expression for this in terms of the bandwidths by
applying Taylor series expansions to both the functional g and
the densities.

The Taylor series expansion of g( X(𝐗) Y(𝐘)ν_1ν_2/ Z(𝐗,𝐘)ν_3)
around f_X(𝐗)f_Y(𝐘)ν_1ν_2 and f_XY(𝐗,𝐘)ν_3
gives an expansion with terms of the form of

Define g(t_1/t_2)=g(t_1,t_2). Define
the following terms:


    Z^q(𝐙)     =ν_3^q( Z(𝐙)-f_XY(𝐗,𝐘))^q,
     XY^q(𝐙)     =(ν_1ν_2)^q( X(𝐗) Y(𝐘)-f_X(𝐗)f_Y(𝐘))^q.

Then the Taylor series expansion of g( X(𝐗) Y(𝐘)ν_1ν_2/ Z(𝐗,𝐘)ν_3)
around f_X(𝐗)f_Y(𝐘)ν_1ν_2 and f_XY(𝐗,𝐘)ν_3
gives

    g( X(𝐗) Y(𝐘)ν_1ν_2/ Z(𝐗,𝐘)ν_3)=∑_i=0^∞∑_j=0^∞(.∂^i+jg(t_1,t_2)/∂ t_1^i∂ t_2^j|_[ t_1=f_X(𝐗)f_Y(𝐘)ν_1ν_2;       t_2=f_XY(𝐗,𝐘)ν_3 ]) Z^j(𝐙) XY^i(𝐙)/i!j!.

To simplify this, we focus on the terms in (<ref>). Note
that if ν_i=1, then the terms in (<ref>) are unaffected.
For other values, ν_i^j decreases to zero as j→∞
since 0<ν_i<1. 

By the binomial theorem, we obtain 

    Z^q(𝐙)     =ν_3^q∑_k=0^q( Z(𝐙))^k(f_XY(𝐗,𝐘))^q-k(-1)^q-k,
     XY^q(𝐙)     =(ν_1ν_2)^q∑_k=0^q( X(𝐗) Y(𝐘))^k(f_X(𝐗)f_Y(𝐘))^q-k(-1)^q-k.

To derive the bias, we will take the conditional expectation given
 of these terms.

Since we are not doing explicit boundary correction, we need to consider
separately the cases when 𝐙 is in the interior of the
support 𝒮_X×𝒮_Y and when 𝐙
is close to the boundary of the support. For precise definitions,
a point Z=(X,Y)∈𝒮_X×𝒮_Y is in the
interior of 𝒮_X×𝒮_Y if for all Z^'∉𝒮_X×𝒮_Y,
K_X(X-X^'/h_X)K_Y(Y-Y^'/h_Y)=0,
and a point Z∈𝒮_X×𝒮_Y is near the
boundary of the support if it is not in the interior. 



 §.§ Interior Points


We first consider the case where 𝐙=(𝐗,𝐘)
is drawn from f_XY in the interior of 𝒮_X×𝒮_Y.
It can be shown (see <cit.>) by Taylor series expansions
of the probability densities that 


    X[ X(𝐗)]     =     f_X(𝐗)+∑_j=1^⌊ s/2⌋c_X,j(𝐗)h_X^2j+O(h_X^s),
     Y[ Y(𝐘)]     =     f_Y(𝐘)+∑_j=1^⌊ s/2⌋c_Y,j(𝐘)h_Y^2j+O(h_Y^s),
    [ Z(𝐙)]     =     f_XY(𝐗,𝐘)+∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j(𝐗,𝐘)h_X^2ih_Y^2j+O(h_X^s+h_Y^s).

The constants in the above expressions are independent of the bandwidths
h_X and h_Y and only depend on the densities and their derivatives.

We will also need an expression for Z[ X(𝐗) Y(𝐘)]:

    Z[ X(𝐗) Y(𝐘)]     =[1/N^2h_X^d_Xh_Y^d_Y∑_i=1^N∑_j=1^NK_X(-_i/h_X)K_Y(-_j/h_Y)]
        = Z[1/N^2h_X^d_Xh_Y^d_Y∑_i=1^NK_X(-_i/h_X)K_Y(-_i/h_Y)]
        + Z[1/N^2h_X^d_Xh_Y^d_Y∑_i,j=1
    
    i≠ j
    ^NK_X(-_i/h_X)K_Y(-_j/h_Y)]
        =1/N Z[ Z(𝐙)]+N^2-N/N^2 X[ X(𝐗)] Y[ Y(𝐘)]
        =N^2-N/N^2f_X(𝐗)f_Y(𝐘)+∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_X,Y,i,j(𝐗,𝐘)h_X^2ih_Y^2j+O(h_X^s+h_Y^s+1/N),

where we have used the fact that _i and _j are independent
when i≠ j.

We will also need expressions for the conditional expectation of powers
of the KDEs to simplify (<ref>). Consider first ( Z(𝐙))^2.
Note that 

    ( Z(𝐙))^2=1/N^2h_X^2d_Xh_Y^2d_Y∑_i=1^N∑_j=1^NK_X(-_i/h_X)K_Y(-_i/h_Y)K_X(-_j/h_X)K_Y(-_j/h_Y).

The above double sum can be split into two cases: i=j (N terms)
and i≠ j (N^2-N terms). When i=j, we have by Taylor
series expansions of the joint density f_XY:

    1/Nh_X^2d_Xh_Y^2d_Y Z[K_X^2(-_i/h_X)K_Y^2(-_i/h_Y)]     =1/Nh_X^d_Xh_Y^d_Y∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,1(𝐗,𝐘)h_X^2ih_Y^2j+O(h_X^s+h_Y^s).


When i j,recall that _i is independent of _j.
Therefore, we obtain for these terms

    N^2-N/N^2h_X^2d_Xh_Y^2d_Y[K_X(-_i/h_X)K_Y(-_i/h_Y)K_X(-_j/h_X)K_Y(-_j/h_Y)]
    
    =N^2-N/N^2 Z[ Z(𝐙)]^2.

Combining these results gives 

    Z[( Z(𝐙))^2]     =N^2-N/N^2f_XY()^2+∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,2(𝐗,𝐘)h_X^2ih_Y^2j
        +∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,3(,)/Nh_X^d_Xh_Y^d_Yh_X^2ih_Y^2j+O(h_X^s+h_Y^s).


For the cross term (k=1) in (<ref>) when q=2, we
obtain 

    Z[ Z(𝐙)]f_XY()=f_XY()^2+f_XY()∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j(𝐗,𝐘)h_X^2ih_Y^2j+O(h_X^s+h_Y^s).

Combining these results gives 

    Z[ Z^2(𝐙)]     =∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,4(𝐗,𝐘)h_X^2ih_Y^2j+∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,5(,)/Nh_X^d_Xh_Y^d_Yh_X^2ih_Y^2j
        -f_XY()^2/N+O(h_X^s+h_Y^s).


By following similar procedures, it can be shown that for q≥2

    Z[ Z^q(𝐙)]     =∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,q,1(𝐗,𝐘)h_X^2ih_Y^2j+∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,q,2(,)/(Nh_X^d_Xh_Y^d_Y)^q-1h_X^2ih_Y^2j
        +O(1/N^q-1+h_X^s+h_Y^s).

In particular, the terms related to f_XY()^q all combine
to be O(1/N). In the q=2 example above, we end up with 

    f_XY()^2-2f_XY()^2+N^2-N/N^2f_XY()^2=-f_XY()^2/N.

As another example, for q=3, we end up with 

    (N(N-1)(N-2)/N^3-3N^2-N/N^2+3-1)f_XY()^3    =(N^3-2N^2-N^2+2N-N^3/N^3-3N^2-N-N^2/N^2)f_XY()^3
        =(-3N+2+3N/N^2)f_XY()^3
        =O(1/N^2).

A similar pattern holds for higher values of q. 

Then by following a similar process, we obtain for q≥2

    Z[ XY^q(𝐙)]     =∑_i=0
    
    i+j≠0
    ^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,q,3(𝐗,𝐘)h_X^2ih_Y^2j+∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,q,4(,)/(Nh_X^d_X)^q-1h_X^2ih_Y^2j
        +∑_i=0^⌊ s/2⌋∑_j=0^⌊ s/2⌋c_XY,i,j,q,5(,)/(Nh_Y^d_Y)^q-1h_X^2ih_Y^2j+O(1/N^q-1+h_X^s+h_Y^s).





 §.§ Points Near the Boundary


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

    X[ X(𝐗)]-f_X()     =1/h_X^d_X∫_V:V∈𝒮_XK_X(𝐗-V/h_X)f_X(V)dV-f_X(𝐗)
        =[1/h_X^d_X∫_V:K_X(𝐗-V/h_X)>0K_X(𝐗-V/h_X)f_X(V)dV-f_X(𝐗)]
        -[1/h_X^d_X∫_V:V∉𝒮_XK_X(𝐗-V/h_X)f_X(V)dV]
        =T_1,X(𝐗)-T_2,X(𝐗).


Note that we are technically evaluating the density f_X at points
outside of the support in T_1,X(). However, to obtain an expression
for this integral, we take a Taylor series expansion of f_X at
the point  which is inside the support. Thus the exact manner
in which we define the extension of f_X 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 T_1,X() gives an expression
similar to that of the interior point case in (<ref>).

For the T_2,X(𝐗) term, we can use multi-index notation
on the expansion of f_X to show that 

    T_2,X(𝐗)     =[1/h_X^d_X∫_V:V∉𝒮_XK_X(𝐗-V/h_X)f_X(V)dV]
        =∫_u:h_Xu+𝐗∉𝒮_X,K_X(u)>0K_X(u)f_X(𝐗+h_Xu)du
        =∑_|α|≤ rh_X^|α|/α!∫_u:h_Xu+𝐗∉𝒮_X,K_X(u)>0K_X(u)D^αf_X(𝐗)u^αdu+o(h_X^r).

Then since the |α|th derivative of f_X is r-|α|
times differentiable, we apply the condition in assumption 𝒜.5
to obtain 

    [T_2,X(𝐗)]=∑_i=1^re_ih_X^i+o(h_X^r).

Similar expressions can be obtained for Y, Z, and the
product Y X.

The above results considers Z^q(𝐙) and XY^q(𝐙)
for q=1. We now consider when q≥2. We follow a similar procedure
where we extend the density beyond the support, but only evaluate
the densities and their derivatives at points within the support.
Thus by the binomial theorem, we can write

    Z[ Z^q(𝐙)]     =ν_3^q∑_k=0^q Z[( Z(𝐙))^k](f_XY(𝐗,𝐘))^q-k(-1)^q-k
        =ν_3^q∑_k=0^q Z[( Z(𝐙))^k]_extended(f_XY(𝐗,𝐘))^q-k(-1)^q-k
        -ν_3^q∑_k=1^q Z[( Z(𝐙))^k]_outside(f_XY(𝐗,𝐘))^q-k(-1)^q-k
        =T_1,q,Z()-T_2,q,Z().

As before, T_1,q,Z() corresponds to the case where we have
extended the density beyond the support and results in terms of the
form in (<ref>). T_2,q,Z() corresponds to the case where
we integrate outside of the boundary. The additional powers applied
to the KDE simply result in terms with the kernel raised to a power
or Z[ Z(𝐙)] raised to a power. By applying
assumption 𝒜.5, we obtain [T_2,q,X(𝐗)]=∑_i=1^r∑_j=1^re_q,ih_X^ih_Y^j+o(h_X^r+h_Y^r).
Similar results are obtained for XY^q(𝐙).

Combining the results for the interior points and points near the
boundary completes the proof.




§ PROOF OF THEOREM <REF> (CONTINUOUS VARIANCE)


Here we prove Theorem <ref>. The
proof uses the Efron-Stein inequality <cit.>: 
 Let 𝐗_1,…,𝐗_n,𝐗_1^',…,𝐗_n^'
be independent random variables on the space 𝒮. Then
if f:𝒮×…×𝒮→ℝ,
we have that 

    [f(𝐗_1,…,𝐗_n)]≤1/2∑_i=1^n[(f(𝐗_1,…,𝐗_n)-f(𝐗_1,…,𝐗_i^',…,𝐗_n))^2].


In this case we consider the samples {𝐙_1,…,𝐙_N}
and {𝐙_1^',_2…,𝐙_N}
and the respective estimators  and ^'. By the triangle
inequality, 

    |-^'|    ≤   1/N|g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^')ν_3)|
           +1/N∑_j=2^N_2|g( X(𝐗_j) Y(𝐘_j)ν_1ν_2/ Z(𝐗_j,𝐘_j)ν_3)-g( X^'(𝐗_j) Y^'(𝐘_1)ν_1ν_2/ Z^'(𝐗_1,𝐘_1)ν_3)|.

By the Lipschitz condition on g, the first term in (<ref>)
can be decomposed into terms of the form of 

    ν_3| Z(𝐙_1)- Z(𝐙_1^')|,


    ν_1ν_2| X(𝐗_1) Y(𝐘_1)- X(𝐗_1^') Y(𝐘_1^')|.


For the Z(𝐙_1) term, we first apply Jensen's inequality:


    [| Z(𝐙_1)- Z(𝐙_1^')|^2]     =[1/N^2h_X^2d_Xh_Y^2d_Y(∑_i=1^N(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)-(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y)))^2]
       ≤1/Nh_X^2d_Xh_Y^2d_Y∑_i=1^N[(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)-(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y))^2]

By making the substitutions 𝐮_i=_1-_i/h_X,
𝐯_i=_1-_i/h_Y, 𝐮_i^'=_1-_i/h_X,
and 𝐯_i^'=_1^'-_i/h_Y in the expectation,
we obtain

    [1/h_X^2d_Xh_Y^2d_Y(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)-(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y))^2]


    =1/h_X^2d_Xh_Y^2d_Y∫(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)-(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y))^2f_Z(_i)f_Z(_1^')f_Z(_1)d_id_1d_1^'
       ≤2||K_X· K_Y||_∞^2.

This gives 

    [ν_3^2| Z(𝐙_1)- Z(𝐙_1^')|^2]≤2||K_X· K_Y||_∞^2,

where we have used the fact that ν_3≤1.

For the product of the marginal KDEs, we have that 

    X(𝐗_1) Y(𝐘_1)     =    1/N^2h_X^d_Xh_Y^d_Y∑_i=2^N∑_j=2^NK_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)
        =    1/N Z(𝐙_1)+1/N^2h_X^d_Xh_Y^d_Y∑_i≠ jK_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y).

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

    [ν_1^2ν_2^2| X(𝐗_1) Y(𝐘_1)- X(𝐗_1^') Y(𝐘_1^')|^2]    ≤   [2/N^2| Z(𝐙_1)- Z(𝐙_1^')|^2]
           +2(M-1)/N^3h_X^2d_Xh_Y^2d_Y×
          ∑_i≠ j[(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)..
          ..-K_X(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y))^2]
       ≤   4+2(N-1)^2/N^2||K_X· K_Y||^2.


For the second term in (<ref>), it can be shown that
(see <cit.>)

    [ν_3^2| Z(𝐙_i)- Z^'(𝐙_i)|^2]     =    ν_3^2/N^2h_X^2d_Xh_Y^2d_Y[(K_X(𝐗_1-𝐗_i/h_X)K_Y(𝐘_1-𝐘_j/h_Y)..
          ..-K_X(𝐗_1^'-𝐗_i/h_X)K_Y(𝐘_1^'-𝐘_j/h_Y))^2]
       ≤   2||K_X· K_Y||_∞^2/N^2.

By a similar approach, 

    X(𝐗_i) Y(𝐘_i)- X^'(𝐗_i) Y^'(𝐘_i)


    =     Z(𝐙_i)- Z^'(𝐙_i)+1/M^2h_X^d_Xh_Y^d_Y(∑_n=2
    
    n≠ i
    K_Y(𝐘_i-𝐘_n/h_Y)(K_X(𝐗_i-𝐗_1/h_X)-K_X(𝐗_i-𝐗_1^'/h_X)).
          .+∑_n=2
    
    n≠ i
    K_X(𝐗_i-𝐗_n/h_X)(K_Y(𝐘_i-𝐘_1/h_Y)-K_Y(𝐘_i-𝐘_1^'/h_Y))),


    [ν_1^2ν_2^2| X(𝐗_i) Y(𝐘_i)- X^'(𝐗_i) Y^'(𝐘_i)|^2]≤6||K_X· K_Y||_∞^2(1/N^2+(N-2)^2/N^4)


We can then apply the Cauchy Schwarz inequality to bound the square
of the second term in (<ref>) to get 

    [(∑_j=2^N_2|g( X(𝐗_1) Y(𝐘_1)/ Z(𝐗_1,𝐘_1))-g( X^'(𝐗_1) Y^'(𝐘_1)/ Z^'(𝐗_1,𝐘_1))|)^2]≤14C_g^2||K_X· K_Y||_∞^2.

Applying Jensen's inequality in conjunction with these results gives

    [|-^'|^2]≤44C_g^2||K_X· K_Y||_∞^2/N^2.

Applying the Efron-Stein inequality finishes the proof.



§ PROOF OF MINIMAX RATES (THEOREM <REF>)


Here we present a proof of the minimax lower
bound on MI estimation convergence rates given in Theorem <ref>.
We follow a similar approach to that given in <cit.>,
which uses the standard approach of Le Cam's method <cit.>.
However, the minimax theory previously derived in these references
are not directly applicable to the MI estimation setting. On the one
hand, MI estimation could be considered to be similar to the entropy
estimation problem as MI is technically a functional of just the joint
distribution as the marginals are derived from the joint. However,
this viewpoint results in a very complicated functional of the joint
distribution, and it is not at all obvious if the theory derived in
<cit.> is applicable. On the other hand, MI estimation
could be viewed as a divergence estimation problem between the joint
and the product of marginals. However, perturbing the joint density,
as is done in Le Cam's method, should result in a perturbation of
the marginals as well. This is not accounted for in existing approaches
for the divergence estimation problem <cit.>.
Therefore, we tailor Le Cam's method directly to the MI estimation
case to avoid these problems.

Le Cam's method reduces the estimation problem to a testing problem
which can then be used to characterize the minimax rate. We will need
the squared Hellinger distance between two densities p and q,
which is defined as 

    H^2(P,Q)=∫(√(p(x,y))-√(q(x,y)))^2dxdy.


For this proof, we present a shift in notation for clarity. We can
write 

    I(;)=I(f_XY,f_Xf_Y).

We now present Le Cam's method adapted to the continuous MI estimation
setting:

Let I be a functional defined on some subset
of a parameter space Θ×Θ which contains (f_XY,f_Xf_Y)
and (f_XY,λ,f_Xf_Y) for all λ in
some index set Λ. Denote the distributions of f_XY,f_Xf_Y,
and f_XY,λ as F_XY, F_XF_Y, and F_XY,λ,
respectively. Define F̅_XY^N=1/|Λ|∑_λ∈ΛF_XY,λ^N.
Consider the following two conditions for γ<2 and β>0:


    (i)  H^2(F̅_XY^N×(F_XF_Y)^N,F_XY^N×(F_XF_Y)^N)≤γ<2,


    (ii)  I(f_XY,f_Xf_Y)≥2β+I(f_XY,λ,f_Xf_Y) ∀λ∈Λ


    OR I(f_XY,λ,f_Xf_Y)≥2β+I(f_XY,f_Xf_Y) ∀λ∈Λ.

 Then 

    inf_Ĝ_Nsup_f_XY∈Θ[|Ĝ_N-I(f_XY,f_Xf_Y)|>β]≥1/2(1-√(γ(1-γ/4))).


Lemma <ref> is nearly identical to Lemma 7 in Krishnamurthy
et al. <cit.> which is itself a modification
of Theorem 2.2 in Tsybakov <cit.>. The
main change is the addition of the second option to condition (ii).
Thus the only required addition to the proof is to check that the
second option is indeed sufficient. Since we consider the probability
of the absolute value of Ĝ_N-I(f_XY,f_Xf_Y)
being greater than β, it is clear that following the original
proof with the second option will give the same result as the first
option in condition (ii).

Note that we are only required to perturb the joint density f_XY
to derive the minimax rate. In general, perturbing the joint density
will result in perturbed marginal densities as well. However, the
following Lemma shows that it is possible to construct perturbation
functions that do not affect the marginal densities. 

Let R_1,…,R_ℓ be a partition of
[0,1]^d_X+d_Y, each being cubes with side length ℓ^-1/(d_X+d_Y).
Then there exists functions u_1,…,u_ℓ such that, 

    supp(u_j)⊂{.z|B(z,ϵ)⊂ R_j} ,


    ∫ u_j^2(x,y)dxdy∈Θ(ℓ^-1),


    ∫ u_j(x,y)dx=∫ u_j(x,y)dy=0,


    ∫ g'(f_X(x)f_Y(y)/f_XY(x,y))f_X(x)f_Y(y)u_j(x,y)/f_XY(x,y)dxdy=0,


    ‖ D^ru_j‖ _∞≤ℓ^r/(d_X+d_Y) ∀ r s.t. ∑_jr_j≤ s+1,

where B(z,ϵ) denotes an L_2 ball around z=[[ x; y ]] with radius ϵ∈(0,1).


Let ϵ>0. We can construct two orthonormal systems of q>3
functions. Construct the first system on [0,1]^d_X such that
ϕ_X,1=1, supp(ϕ_X,j)⊂[ϵ,1-ϵ]^d_X,
and ‖ D^rϕ_X,j‖ _∞≤ J_X<∞
for all j. The second system is constructed on [0,1]^d_Ysuch
that ϕ_Y,1=1, supp(ϕ_Y,j)⊂[ϵ,1-ϵ]^d_Y,
and ‖ D^rϕ_Y,j‖ _∞≤ J_Y<∞
for all j. We can then construct a combined orthonormal system
ϕ_i,j=ϕ_X,iϕ_Y,j which has q^2 functions. It
is clear that this is an orthonormal system since ∫ϕ_i,j(x,y)ϕ_m,n(x,y)dxdy=1
when i=m and j=n and zero otherwise. Also, ϕ_1,1=1,
supp(ϕ_i,j)⊂[ϵ,1-ϵ]^d_X+d_Y,
and there exists some J<∞ such that ‖ D^rϕ_i,j‖ _∞≤ J
for all i,j.

Now for any given function f∈ L_2([0,1]^d_X+d_Y),
we can find a unit-normed function v∈span({ϕ_i,j})
such that v⊥ϕ_1,i for all i, v⊥ϕ_i,1 for
all i, and v⊥ f. We can write v=∑_i,j=1^qa_ib_jϕ_i,j.
Then D^rv=∑_i,j=1^qa_ib_jD^rϕ_i,j which implies
that 

    ‖ D^rv‖ _∞≤ J∑_i,j|a_ib_j|≤ Jq,

where the last inequality comes from the fact that v is unit-normed.
Now define ν=1/Jqv. Then clearly ∫ν^2(x,y)dxdy
is upper and lower bounded and we have that ‖ D^rν‖ _∞≤1.

We can now construct the functions u_j. First map R_j to
[0,1]^d_X+d_Yby scaling it appropriately. Then set 

    u_j(x,y)=ν(ℓ^1/(d_X+d_Y)([[ x; y ]]-𝐣)),

where 𝐣 is the point in R_j that is mapped to 0
after scaling. This maps [0,1]^d_X+d_Y back to R_j while
inheriting the properties derived from the construction of ν.
Now let f be g'(f_X(x)f_Y(y)/f_XY(x,y))f_X(x)f_Y(y)/f_XY(x,y)
constrained to R_j and scaled to fit [0,1]^d_X+d_Y.
Conditions 1, 3, and 4 above are then fulfilled by construction. Also
∫_R_ju_j^2(x,y)dxdy=1/ℓν^2(x,y)dxdy∈Θ(ℓ^-1)
which is condition 2 above. It's also clear that ‖ D^ru_j‖ _∞≤ℓ^r/(d_X+d_Y)
which is condition 5, completing the proof.

We now can prove Theorem <ref>. We will construct the
conditions necessary to apply Lemma <ref>. Apply Lemma <ref>
to obtain an index set Λ̃={-1,1}^ℓ and functions
u_1,…,u_ℓ. Define the following set of perturbed functions
around f_XY:

    Λ={ f_XY,λ=f_XY+H_1∑_j=1^ℓλ_ju_j|λ_j∈Λ̃} .

This will form our set of alternatives. Due to the fact that ∫ u_j(x,y)dx=∫ u_j(x,y)dy=0,
we have that 

    ∫ f_XY,λ(x,y)dx=f_Y(y),


    ∫ f_XY,λ(x,y)dy=f_X(x).

That is, the perturbations on f_XY are chosen so that the resulting
marginal distributions are unperturbed. 

The perturbation functions u_j in Lemma <ref>
are restricted to the small R_j bins and thus violate the Hölder
class assumption. However, by scaling H_1 appropriately, we can
ensure that f_XY,λ∈Σ(s,H). We show this by following
the same argument as in Krishnamurthy et al. <cit.>,
which we repeat here for completeness. Define u_λ=H_1∑_j=1^ℓλ_ju_j.
We will first show that u_λ is Hölder smooth. Then
f_XY,λ is Hölder smooth by the triangle inequality.
For u_λ, fix two points v,z∈ℝ^d_X+d_Y
and fix r with ∑_jr_j=s. Define z_1 as the boundary
point of the R_j bin containing z along the line between z
and v and define v_1 similarly as the boundary point for the
bin containing v along the same line. We then have the following:

    |D^ru_λ(z)-D^ru_λ(v)|    ≤|D^ru_λ(z)-D^ru_λ(z_1)|+|D^ru_λ(z_1)-D^ru_λ(v_1)|+|D^ru_λ(v_1)-D^ru_λ(v)|
        =|D^ru_λ(z)-D^ru_λ(z_1)|+|D^ru_λ(v_1)-D^ru_λ(v)|
        =∫_γ(z,z_1)∇ D^ru_λ(t)dt+∫_γ(v,v_1)∇ D^ru_λ(t)dt
       ≤ H_1‖ D^r+1u_j‖ _∞(‖ z-z_1‖ _2+‖ v-v_1‖ _2)
       ≤ H_1ℓ^r+1/d_X+d_Y(‖ z-z_1‖ _2^1-(s-r)‖ z-z_1‖ _2^s-r+‖ v-v_1‖ _2^1-(s-r)‖ v-v_1‖ _2^s-r)
       ≤ H_1ℓ^r+1/d_X+d_Y√(d_X+d_Y)ℓ^-1-(s-r)/d_X+d_Y(‖ z-z_1‖ _2^s-r+‖ v-v_1‖ _2^s-r)
       ≤ H_1ℓ^s/d_X+d_Y√(d_X+d_Y)‖ z-v‖ _2^s-r.


The first line is an application of the triangle inequality. The second
line follows from the fact that u_λ and all of its derivatives
are zero on the boundaries of the cubes R_j as u_j is not
supported in the band around the border of R_j. The third line
follows from the fundamental theorem of calculus where γ(z,z_1)
is the path between z and z_1. The fourth line is an application
of Hölder's inequality where we replace each derivative with
its supremum, leaving just the path integral which simplifies to the
length of the path. The fifth line follows from the assumption that
‖ D^ru_j‖ _∞≤ℓ^r/(d_X+d_Y)
when ∑_jr_j≤ s+1. For the sixth line, since z and
z_1 are in the same bin, then ‖ z-z_1‖ _2≤√(d_X+d_Y)ℓ^-1/(d_X+d_Y)
as there are ℓ boxes with side length ℓ^-1/(d_X+d_Y).
Finally, the last line follows since z_1 and v_1 are on
the line segment between z and v.

This indicates that u_λ, and therefore f_XY,λ,
is guaranteed to be Hölder smooth if H_1ℓ^s/d_X+d_Y√(d_X+d_Y)≤ H.
Thus we require that H_1=O(ℓ^-s/d_X+d_Y).
We will set ℓ later on.

Note that for any f_XY,λ∈Λ, by a second order Taylor
series approximation in the first argument we have 

    I(f_XY,λ,f_Xf_Y)     =I(f_XY,f_Xf_Y)-∫ g'(f_X(x)f_Y(y)/f_XY(x,y))f_X(x)f_Y(y)u_λ(x,y)/f_XY(x,y)dxdy
        +1/2∫ g”(f_X(x)f_Y(y)/f_XY^*(x,y))f_X^2(x)f_Y^2(y)u_λ^2(x,y)/(f_XY^*(x,y))^2dxdy,

where f_XY^* is the function from Taylor's remainder theorem.
By construction (see Lemma <ref>), the first order
term vanishes. For the second order term, note that f_XY^*
lies on the line segment between f_XY and f_XY,λ and
is therefore upper and lower bounded. Similarly the density f_XY,λ
will be upper and lower bounded for N sufficiently large as f_XY,λ∈[f_XY-H_1,f_XY+H_1]
due to the fact that ‖ D_0u_j‖ _∞=‖ u_j‖ _∞≤1,
and H_1 will be chosen to decrease as N increases. Assume
without loss of generality that given ϵ>0, g”(ϵ)>0.
Thus there exists a constant c_0 such that 

    1/2∫ g”(f_X(x)f_Y(y)/f_XY^*(x,y))f_X^2(x)f_Y^2(y)u_λ^2(x,y)/(f_XY^*(x,y))^2dxdy≥ c_0H_1^2∑_j=1^ℓ‖ u_j‖ _2^2≥ c_1H_1^2,

where we have used the facts that the u_j functions are orthogonal
to each other and ‖ u_j‖ _2^2∈Θ(ℓ^-1).
Therefore,

    I(f_XY,λ,f_Xf_Y)-I(f_XY,f_Xf_Y)≥ c_1H_1^2,

providing us with the necessary separation of 2β where β=c_1H_1^2/2.
Note that if g”(ϵ)<0, we can simply consider I(f_XY,f_Xf_Y)-I(f_XY,λ,f_Xf_Y)
instead.

We now focus on bounding the squared Hellinger distance H^2(F̅_XY^N×(F_XF_Y)^N,F_XY^N×(F_XF_Y)^N)
where F̅_XY^N=1/|Λ|∑_λ∈ΛF_XY,λ^N.
The Hellinger distance decomposes across product measures resulting
in:

    H^2(F̅_XY^N×(F_XF_Y)^N,F_XY^N×(F_XF_Y)^N)     =2(1-(1-H^2(F̅_XY^N,F_XY^N)/2)(1-H^2((F_XF_Y)^N,(F_XF_Y)^N)))
        =H^2(F̅_XY^N,F_XY^N).

To bound this, we will use the following result from Birge and Massart <cit.>:

Consider a set of densities p and p_λ=p(1+∑_jλ v_j)
for λ∈Λ={-1,1}^ℓ. Suppose that (i) ‖ v_j‖ _∞≤1,
(ii) ‖ 1_{R_j^C}v_j‖ _1=0, (iii)
∫ pv_j=0, and (iv) ∫ pv_j^2=α_j>0 all hold.
Define P̅^N=1/|Λ|∑_λ∈ΛP_λ^N.
Then 

    H^2(P̅^N,P^N)≤N^2/3∑_j=1^ℓα_j^2.


To apply Lemma <ref>, define v_j(z)=H_1u_j(z)/f_XY(z).
Then f_XY,λ=f_XY(1+∑_jλ_jv_j).
Requirements (i)-(iii) are immediately satisfied based on the properties
of u_j (see Lemma <ref>). Furthermore, 

    α_j=∫ v_j^2f_XY=H_1^2∫ u_j^2/f_XY≤H_1^2C/ℓ,

for some constant C. Therefore 

    H^2(F̅_XY^N,F_XY^N)≤N^2/3∑_j=1^ℓα_j^2≤N^2H_1^4C^2/ℓ∈Θ(N^2ℓ^-4s+d_X+d_Y/d_X+d_Y).

Set ℓ=N^2(d_X+d_Y)/4s+d_X+d_Y resulting
in H_1=N^-2s/4s+d_X+d_Y. Then the Hellinger distance
is bounded by a constant. Additionally, the error is larger than β∈Θ(N^-4s/4s+d_X+d_Y)
allowing us to apply Lemma <ref> when s<(d_X+d_Y)/4.
Markov's inequality then finishes the proof.

For s>(d_X+d_Y)/4,we get a lower bound of O(N^-1)
which is the parametric rate. In general, we cannot do any better
than this <cit.>
thus establishing the lower bound in this regime. In particular, Krishnamurthy
et al. <cit.> use a contradiction approach
to establish this for divergence estimation which can be extended
to the MI estimation problem.



§ THEORY FOR MIXED RANDOM VARIABLES


Here we provide proofs of the theory that extends the MI estimators
for the continuous case to the mixed case.



 §.§ Proof of Lemma <ref>




For (<ref>), note that_xy is a binomial random
variable with parameter f_X_DY_D(x,y), N trials,
and mean Nf_X_DY_D(x,y). Thus (<ref>)
is the (potentially) fractional moment of a binomial random variable.
By the generalized binomial theorem, we have that

    𝐍_xy^α    =    (𝐍_xy-Nf_X_DY_D(x,y)+Nf_X_DY_D(x,y))^α
        =    ∑_i=0^∞([ α; i ])(Nf_X_DY_D(x,y))^α-i(𝐍_xy-Nf_X_DY_D(x,y))^i,
    [𝐍_xy^α]     =    ∑_i=0^∞([ α; i ])(Nf_X_DY_D(x,y))^α-i[(𝐍_xy-Nf_X_DY_D(x,y))^i].

From <cit.>, the i-th central moment of 𝐍_xy
has the form of 

    [(𝐍_xy-Nf_X_DY_D(x,y))^i]=∑_n=0^⌊ i/2⌋c_n,i(f_X_DY_D(x,y))N^n.

Combining this with (<ref>) gives 

    [_xy^α]     =∑_i=0^∞∑_n=0^⌊ i/2⌋([ α; i ])(f_X_DY_D(x,y))^α-ic_n,i(f_X_DY_D(x,y))N^α-i+n
        =(Nf_X_DY_D(x,y))^α+O(N^α-1).


For (<ref>), we apply a Taylor series expansion
to obtain

    _xy^λ_x^β_y^γ    =N^λ+β+γp^λp_x^βp_y^γ+(_xy-Np)p^λ-1(N^λ+β+γ-1p_x^βp_y^γ+N^λ+β+γ-2(p_x^β-1p_y^γ(_x-Np_x)+p_x^βp_y^γ-1(_y-Np_y)))
        +N^λ+β+γ-1p^λ(p_x^β-1p_y^γ(_x-Np_x)+p_x^βp_y^γ-1(_y-Np_y))+O(N^λ+β+γ-2((_x-Np_x)(_y-Np_y))),

where we set p=f_X_DY_D(x,y), p_x=f_X_D(x),
and p_y=f_Y_D(y) for notational convenience. By taking the
expected value with respect to _x, _y, and _xy,
we obtain 

    [_xy^λ_x^β_y^γ]     =N^λ+β+γp^λp_x^βp_y^γ+N^λ+β+γ-2p^λ-1(p_x^β-1p_y^γ(_xy,_x)+p_x^βp_y^γ-1(_xy,_y))
        +O(N^β+γ-1(_x,_y))
        =N^λ+β+γp^λp_x^βp_y^γ+O(N^λ+β+γ-1),

where the last step follows from the Cauchy-Schwarz inequality and
the variance of a binomial random variable.



 §.§ Proof of Theorem <ref> (Bias)




For notational ease, let 

    𝒯(,)=f_X_C|X_D(_C|_D)f_Y_C|Y_D(_C|_D)/f_X_CY_C|X_DY_D(_C,_C|_D,_D).

We have that 

    [h_X_C|X_D,h_Y_C|Y_D]     =[h_X_C|X_D,h_Y_C|Y_D]-I(;)
        =[∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/Nh_X_C|x,h_Y_C|y-g(𝒯(,)×f_X_D(_D)f_Y_D(_D)/f_X_DY_D(_D,_D))]
        =[∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/N(h_X_C|x,h_Y_C|y-g(𝒯(,)×_x_y/N_xy))]
        +[∑_x∈𝒮_X_D,y∈𝒮_Y_D(_xy/Ng(𝒯(,)×_x_y/N_xy)-f_X_DY_D(x,y)g(𝒯(,)×f_X_D(x)f_Y_D(y)/f_X_DY_D(x,y)))].

We consider the second term in (<ref>) first. A Taylor
series expansion of g(𝒯(,)×_x_y/N_xy)
evaluated at 𝒯(,)×f_X_D(x)f_Y_D(y)/f_X_DY_D(x,y)
gives terms of the form of 

    (f_X_C|X_D(_C|x)f_Y_C|Y_D(_C|y)(_x_y/N^2-f_X_D(x)f_Y_D(y)))^i,


    (f_X_CY_C|X_DY_D(_C,_C|x,y)(_xy/N-f_X_DY_D(x,y)))^i,

where i is a positive integer. For notational ease, set p=f_X_DY_D(x,y).
By applying the binomial theorem and (<ref>), we obtain

    _xy/N(p-_xy/N)^i    =∑_k=0^iikp^i-k(_xy/N)^k+1(-1)^k
    [_xy/N(p-_xy/N)^i]     =p^i+1∑_k=0^iik(-1)^k+O(1/N)
        =O(1/N).


Using a similar approach with (<ref>), it can be
shown that 

    [_xy/N(_x_y/N^2-f_X_D(x)f_Y_D(y))^i]=O(1/N).

Thus the second term in (<ref>) reduces to O(1/N).

By conditioning on _1,D,…,_N,D,_1,D,…,_N,D,
the first term in (<ref>) can be written as 

    [∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/N[.h_X_C|x,h_Y_C|y|_1,D,…,_N,D,_1,D,…,_N,D]].

The conditional bias of h_X_C|x,h_Y_C|y given _1,D,…,_N,D,_1,D,…,_N,D
can be obtained from Theorem <ref> as 

    [.h_X_C|x,h_Y_C|y|_1,D,…,_N,D,_1,D,…,_N,D]     =∑_i,j=0
    
    i+j≠0
    ^rc_10,i,j(_x_y/N^2,_xy/N)^i^j
        +O(^s+^s+1/_xy^d_X^d_Y).

This expression provides the motivation for our choice of  and
. Since ∝_x^-β and ∝_y^-α,
then (<ref>) gives terms with the form of _xy_x^-β i_y^-α j/N
with i+j≥1. From Lemma <ref>, taking
the expected value of these terms gives

    [_xy_x^-β i_y^-α j/N]=N^-β i-α jf_X_DY_D(x,y)(f_X_D(x))^-β i(f_Y_D(y))^-α j+o(1/N).

Similarly, taking the expectation of _xy_x^β d_X_y^α d_Y/N^2
gives O(N^β d_X+α d_Y-1). Note that the
polynomial terms of _x_y/N^2 and _xy/N in the
constants in (<ref>) do not contribute to the
bias rate as the _x_y and _xy terms in the numerator
are cancelled by the N^2 and N terms in the denominator, respectively,
after taking the expectation. Combining all of these results completes
the proof.



 §.§ Proof of Theorem <ref> (Variance)


By the law of total variance, which can
be derived from the Pythagorean theorem, we have 

    [h_X_C|X_D,h_Y_C|Y_D]     =[[.h_X_C|X_D,h_Y_C|Y_D|_1,D,…,_N,D,_1,D,…,_N,D]]
        +[[.h_X_C|X_D,h_Y_C|Y_D|_1,D,…,_N,D,_1,D,…,_N,D]].

Given all of the _i,D and _i,D random variables, the
estimators h_X_C|x,h_Y_C|y are all conditionally independent
since they use different sets of 𝐗_i,C's and _i,C's
for each pair (x,y). Thus from Theorem <ref>, we
get 

    [.h_X_C|X_D,h_Y_C|Y_D|_1,D,…,_N,D,_1,D,…,_N,D]     =O(∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy^2/N^21/_xy)
        =O(∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/N^2).

 Taking the expectation yields O(1/N).

For the second term in (<ref>), we know from (<ref>)
that 

    [.h_X_C|x,h_Y_C|y|_1,D,…,_N,D,_1,D,…,_N,D]     =O(∑_i,j=0
    
    i+j≠0
    ^r_x^-iβ_y^-jα+_x^-sβ+_y^-sα+_x^β d_X_y^α d_Y/_xy)
        =O(f(_x,_y,_xy)).

Let _xy', _x', and _y' be independent and identically
distributed realizations of _xy, _x, and _y,
respectively. Then by the Efron-Stein inequality, 

    [∑_x∈𝒮_X_D,y∈𝒮_Y_D_xy/Nf(_x,_y,_xy)]    ≤1/2N^2∑_x∈𝒮_X_D,y∈𝒮_Y_D[(_xyf(_x,_y,_xy)-_xy'f(_x',_y',_xy'))^2],

where since _x, _y, and _xy are not independent,
we consider the effect of resampling all three simultaneously. Note
that 

    (_xyf(_x,_y,_xy)-_xy'f(_x,'_y',_xy'))^2    =O((∑_i,j=0
    
    i+j≠0
    ^r(_xy_x^-iβ_y^-jα-_xy'(_x')^-iβ(_y')^-jα)..
        +(_xy_x^-sβ-_xy'(_x')^-sβ)+(_xy_y^-sα-_xy'(_y')^-sα)
       ..[ ; ; ; ;  ]+(_x^β d_X_y^α d_Y-(_x')^β d_X(_y')^α d_Y))^2).

By Jensen's inequality, we can consider separately each of the squared
differences in (<ref>). Then since (_xy,_x,_y)
is independent of (_xy',_x',_y') and they
are identically distributed, then the expected squared difference
is proportional to the variance. For example, applying Lemma <ref>
gives

    [(_xy_x^-sβ-_xy'(_x')^-sβ)^2]     =2[_xy_x^-sβ]
        =2N^2-2sβ(f_X_DY_D(x,y))^2(f_X_D(x))^-2sβ-2(N^1-sβf_X_DY_D(x,y)(f_X_D(x))^-sβ)^2
        +O(N^1-2sβ)
        =O(N^1-2sβ).


By a similar procedure, we obtain 

    [(_xy_y^-sα-_xy'(_y')^-sα)^2]     =O(N^1-2sα),
    [(_x^β d_X_y^α d_Y-(_x')^β d_X(_y')^α d_Y)^2]     =O(N^2β d_X+2α d_Y-1),
    [(∑_i,j=0
    
    i+j≠0
    ^r(_xy_x^-iβ_y^-jα-_xy'(_x')^-iβ(_y')^-jα))^2]     =O(N^1-2β)+O(N^1-2α).

Combining these results with (<ref>) and (<ref>)
completes the proof.



§ PROOF OF THEOREM <REF> (CLT)


 We will first find the asymptotic distribution
of 

    √(N)(-[])     =1/√(N)∑_i=1^N(g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-Z_i[g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)])
        +1/√(N)∑_i=1^N(Z_i[g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)]-[g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)]).

By the standard central limit theorem <cit.>,
the second term converges in distribution to a Gaussian random variable
with variance 

    [ Z[g( X(𝐗) Y(𝐘)ν_1ν_2/ Z(𝐗,𝐘)ν_3)]].

All that remains is to show that the first term converges in probability
to zero as Slutsky's theorem <cit.> can then be
applied. Denote this first term as _N and note that [_N]=0. 

We will use Chebyshev's inequality combined with the Efron-Stein inequality
to bound the variance of _N. Consider the samples {_1,…,_N}
and {_1^',_2,…,_N} and the respective
sequences _N and _N^'. This gives 

    _N-_N^'    =1/√(N)(g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-Z_1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)])
        -1/√(N)(g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^'))-Z_1^'[g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^'))])
        +1/√(N)∑_i=2^N(g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)).


Note that 

    [(g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-Z_1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)])^2]=[__1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)]].

We will use the Efron-Stein inequality to bound __1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)].
We thus need to bound the conditional expectation of the term 

    |g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-g( X^'(𝐗_1) Y^'(𝐘_1)ν_1ν_2/ Z^'(𝐗_1,𝐘_1)ν_3)|^2,

where _i is replaced with _i^' in the KDEs for some
i≠1. Using similar steps as in Appendix <ref>,
we have that 

    [|g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-g( X^'(𝐗_1) Y^'(𝐘_1)ν_1ν_2/ Z^'(𝐗_1,𝐘_1)ν_3)|^2]=O(1/N^2).

Then by the Efron-Stein inequality, __1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)]=O(1/N).
Therefore 

    [1/N(g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-Z_1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)])^2]=O(1/N^2).

A similar result holds for the g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^')ν_3)
term in (<ref>). 

For the third term in (<ref>), 

    [(∑_i=2^N|g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)|)^2]
        =∑_i,j=2^N[|g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)|.
       ×.|g( X(𝐗_j) Y(𝐘_j)ν_1ν_2/ Z(𝐗_j,𝐘_j)ν_3)-g( X^'(𝐗_j) Y^'(𝐘_j)ν_1ν_2/ Z^'(𝐗_j,𝐘_j)ν_3)|]

For the N-1 terms where i=j, we know from Appendix <ref>
that 

    [|g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)|^2]=O(1/N^2).

Thus these terms contribute O(1/N). For the (N-1)^2-(N-1)
terms where i≠ j, we can do multiple substitutions of the form
𝐮_j=_j-_1/h_X resulting in

    [|g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i))-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i))|.
    ×.|g( X(𝐗_j) Y(𝐘_j)ν_1ν_2/ Z(𝐗_j,𝐘_j)ν_3)-g( X^'(𝐗_j) Y^'(𝐘_j)ν_1ν_2/ Z^'(𝐗_j,𝐘_j)ν_3)|]     =O(h_X^2d_Xh_Y^2d_Y/N^2).

Since h_X^d_Xh_Y^d_Y=o(1), 

    [(∑_i=2^N|g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)|)^2]=o(1).

Combining all of these results with Jensen's inequality gives 

    [(_N-_N^')^2]    ≤3/N[(g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)-Z_1[g( X(𝐗_1) Y(𝐘_1)ν_1ν_2/ Z(𝐗_1,𝐘_1)ν_3)])^2]
        +3/N[(g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^')ν_3)-Z_1^'[g( X(𝐗_1^') Y(𝐘_1^')ν_1ν_2/ Z(𝐗_1^',𝐘_1^')ν_3)])^2]
        +3/N[(∑_i=2^N(g( X(𝐗_i) Y(𝐘_i)ν_1ν_2/ Z(𝐗_i,𝐘_i)ν_3)-g( X^'(𝐗_i) Y^'(𝐘_i)ν_1ν_2/ Z^'(𝐗_i,𝐘_i)ν_3)))^2]
        =o(1/N).

Applying the Efron-Stein inequality gives that [_N]=o(1).
Then by ChebyShev's inequality, _N converges to zero in probability.
This completes the proof for the plug-in estimator.

For the weighted ensemble estimator, we present a more general result
where we have different parameters l_X∈ℒ_X and
l_Y∈ℒ_Y for  and , respectively. We
can then write 

    √(N)( w-[ w])     =1/√(N)∑_i=1^N∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)(g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3).
       .-Z_i[g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3)])
        +1/√(N)∑_i=1^N(Z_i[∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3)].
       .-[∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3)]).

By the central limit theorem, the second term converges in distribution
to a zero-mean Gaussian random variable with variance 

    [Z_i[∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3)]].

From the previous results, the first term converges to zero in probability
as it can be written as 

    ∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)1/√(N)∑_i=1^N(g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3).
    .-Z_i[g( XX(𝐗_i) YY(𝐘_i)ν_1ν_2/ ZZ(𝐗_i,𝐘_i)ν_3)])     =∑_l_X∈ℒ_X,l_Y∈ℒ_Yw(l_X,l_Y)o_P(1)
        =o_P(1),

where o_P(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 <ref> follows a similar
procedure as the extension to the ensemble case.
