---
abstract: 'We derive the mean squared error convergence rates of kernel density-based plug-in estimators of mutual information measures between two multidimensional random variables $\mathbf{X}$ and $\mathbf{Y}$ for two cases: 1) $\X$ and $\Y$ are both continuous; 2) $\X$ is continuous and $\Y$ is discrete. Using the derived rates, we propose an ensemble estimator of these information measures for the second case by taking a weighted sum of the plug-in estimators with varied bandwidths. The resulting ensemble estimator achieves the $1/N$ parametric convergence rate when the conditional densities of the continuous variables are sufficiently smooth. To the best of our knowledge, this is the first nonparametric mutual information estimator known to achieve the parametric convergence rate for this case, which frequently arises in applications (e.g. variable selection in classification). The estimator is simple to implement as it uses the solution to an offline convex optimization problem and simple plug-in estimators. A central limit theorem is also derived for the ensemble estimator. Ensemble estimators that achieve the parametric rate are also derived for the first case ($\X$ and $\Y$ are both continuous) and another case 3) $\X$ and $\Y$ may have any mixture of discrete and continuous components.'
author:
- 'Kevin R. Moon[^1]'
- 'Kumar Sricharan[^2]'
- 'Alfred O. Hero III[^3]'
bibliography:
- 'References.bib'
title: Ensemble Estimation of Mutual Information
---

\#1[\_[\#1,h\_[\#1]{}]{}]{} \#1[\_[\#1]{}]{}

\#1[\_[\#1,h\_[\#1]{}]{}]{} \#1[\_[\#1,h]{}]{} \#1[\_[\#1,h\_[\#1]{}(l)]{}]{} \#1[\_[\#1,h\_[\#1]{}]{}]{} \#1[\_[\#1,h]{}]{}

§ \#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]{}]{}]{}

\#1\#2\#3[\_[\#2]{}\^[(\#3)]{}]{}

Introduction {#sec:intro}
============

Mutual information (MI) estimation has many applications in machine learning including MI has been used in fMRI data processing [@chai2009fMRI], structure learning [@structure2016], independent subspace analysis [@pal2010estimation], forest density estimation [@liu2012exponential], clustering [@lewi2006real], neuron classification [@schneidman2003information], and intrinsically motivated reinforcement learning [@mohamed2015variational; @salge2014changing]. Another particularly common application is feature selection or extraction where features are chosen to maximize the MI between the chosen features $\mathbf{X}$ and the outcome variables $\mathbf{Y}$ [@torkkola2003feature; @vergara2014review; @peng2005feature; @kwak2002input]. In many of these applications, the predictor labels have discrete components (e.g. classification labels) while the input variables have continuous components. To the best of our knowledge, there are currently no nonparametric MI estimators that are known to achieve the parametric mean squared error (MSE) convergence rate $1/N$ when $\X$ and/or $\Y$ contain discrete components. Also, while many nonparametric estimators of MI exist, most can only be applied to specific information measures (e.g. Shannon or Rényi information). In this paper, we provide a framework for nonparametric estimation of a large class of MI measures where we only have available a finite population of i.i.d. samples. **** We separately consider three cases: 1) $\X$ and $\Y$ are both continuous; 2) $\X$ is continuous and $\Y$ is discrete; 3) $\X$ and $\Y$ may have any mixture of discrete and continuous components. We focus primarily on the second case which includes the problem of feature selection in classification. We derive a MI estimator for this case that achieves the parametric MSE rate when the conditional densities of the continuous variables are sufficiently smooth. We also show how these estimators are extended to the first and third cases.

Our estimation method applies to other MI measures in addition to Shannon information, which have been the focus of much interest. The authors of [@torkkola2003feature] defined an information measure based on a quadratic divergence that could be estimated more efficiently than Shannon information. A MI measure based on the Pearson divergence was considered in [@sugiyama2012machine] for computational efficiency and numerical stability. The authors of [@costa2004geodesic] and [@pal2010estimation] used minimal spanning tree generalized nearest-neighbor graph approaches, respectively, to estimate Rényi information.

Related Work
------------

Many estimators for Shannon MI between continuous random variables have been developed. A popular $k$-nn-based estimator was proposed in [@kraskov2004estimating] which is a modification of the entropy estimator derived in [@kozachenko1987sample]. However, these estimators only achieve the parametric convergence rate when the dimension of each of the random variables is less than 3 [@gao2016demystifying]. Similarly, the Rényi information estimator in [@pal2010estimation] does not achieve the parametric rate. Some other estimators are based on maximum likelihood estimation of the likelihood ratio [@suzuki2008approximating] and minimal spanning trees [@khan2007relative].

Recent work has focused on nonparametric divergence estimation for purely continuous random variables. One approach [@krishnamurthy2014divergence; @kandasamy2015nonparametric; @singh2014exponential; @singh2014renyi] uses an optimal kernel density estimator (KDE) to achieve the parametric convergence rate when the densities are at least $d$ [@singh2014exponential; @singh2014renyi] or $d/2$ [@krishnamurthy2014divergence; @kandasamy2015nonparametric] times differentiable where $d$ is the dimension of the data. These optimal KDEs require knowledge of the density support boundary and are difficult to construct near the boundary. Numerical integration may also be required for estimating some divergence functionals under this approach, which can be computationally expensive. In contrast, our approach to MI estimation does not require numerical integration and can be performed without knowledge of the support boundary.

More closely related work [@sricharan2013ensemble; @moon2014isit; @moon2014nips; @moon2016arxiv; @moon2016isit] uses an ensemble approach to estimate entropy or divergence functionals. These works construct an ensemble of simple plug-in estimators by varying the neighborhood size of the density estimators. They then take a weighted average of the estimators where the weights are chosen to decrease the bias with only a small increase in the variance. The parametric rate of convergence is achieved when the densities are either $d$ [@sricharan2013ensemble; @moon2014isit; @moon2014nips] or $(d+1)/2$ [@moon2016arxiv; @moon2016isit] times differentiable. These approaches are simple to implement as they only require simple plug-in estimates and the solution of an offline convex optimization problem. These estimators have also performed well in various applications [@szabo2012distributed; @gliske2015intrinsic; @moon2015Bayes; @moon2015partI; @moon2015partII]

Finally, the authors of [@gao2015efficient] showed that $k$-nn or KDE based approaches underestimate the MI when the MI is large. As MI increases, the dependencies between random variables increase which results in less smooth densities. Thus a common approach to overcome this issue is to require the densities to be smooth [@krishnamurthy2014divergence; @kandasamy2015nonparametric; @singh2014exponential; @singh2014renyi; @sricharan2013ensemble; @moon2014isit; @moon2014nips; @moon2016arxiv; @moon2016isit].

Contributions
-------------

In the context of this related work, we make the following novel contributions in this paper: (1) For continuous random variables (case 1), we extend the asymptotic bias and variance results for divergence estimators [@moon2016isit; @moon2016arxiv] to kernel density plug-in MI estimators without boundary correction [@karunamuni2005boundary] by incorporating machinery to handle the dependence between the product of marginal density estimators (Section \[sec:MI\_est\]), (2) we extend the theory to handle discrete random variables in the mixed cases (cases 2 and 3) by reformulating the densities as a mixture of the conditional density of the continuous variables given the discrete variables (Section \[sec:mixed\]), and (3) we leverage this theory for the mixed cases in conjunction with the generalized theory of ensemble estimators [@moon2016arxiv; @moon2016isit] to derive, to the best of our knowledge, the first non-parametric estimator that achieves a parametric rate of MSE convergence of $O\left(1/N\right)$ for the mixed cases (Section \[sec:mixed\_ensemble\]), where $N$ is the number of samples available from each distribution. We also derive a central limit theorem for the ensemble estimators (Section \[subsec:clt\]). We verify the theory through experiments (Section \[sec:experiments\]).

Continuous Random Variables {#sec:MI_est}
===========================

In this section, we obtain MSE convergence rates of plug-in MI estimators when $\X$ and $\Y$ are continuous (case 1 in Section \[sec:intro\]). This will enable us to derive the MSE convergence rates of plug-in MI estimators when $\X$ is continuous and $\Y$ is discrete and when $\X$ and $\Y$ may have any mixture of continuous and discrete components (respectively, cases 2 and 3 in Section \[sec:intro\]). These rates can then be used to derive ensemble estimators that achieve the parametric MSE rate. Let $f_{X}(x)$, $f_{Y}(y)$, and $f_{XY}(x,y)$ be $d_{X}$, $d_{Y}$, and $d_{X}+d_{Y}=d$-dimensional densities. Let $g(t_{1},t_{2})=g\left(\frac{t_{1}}{t_{2}}\right)$ (e.g. $g(t_{1},t_{2})=\log(t_{1}/t_{2})$ for Shannon information). We define a family of MIs as $$G_{1}(\mathbf{X};\mathbf{Y})=\int g\left(\frac{f_{X}(x)f_{Y}(y)}{f_{XY}(x,y)}\right)f_{XY}(x,y)dxdy.\label{eq:MI}$$

The KDE Plug-in Estimator
-------------------------

When both $\mathbf{X}$ and $\mathbf{Y}$ are continuous with marginal densities $f_{X}$ and $f_{Y}$, the MI functional $G_{1}(\mathbf{X};\mathbf{Y})$ can be estimated using KDEs. Assume that $N$ i.i.d. samples $\left\{ \mathbf{Z}_{1},\dots,\mathbf{Z}_{N}\right\} $ are available from the joint density $f_{XY}$ with $\mathbf{Z}_{i}=\left(\mathbf{X}_{i},\mathbf{Y}_{i}\right)^{T}$. Let $M=N-1$ and let $h_{X}$, $h_{Y}$ be kernel bandwidths. Let $K_{X}(\cdot)$ and $K_{Y}(\cdot)$ be kernel functions with $||K_{X}||_{\infty},\,||K_{Y}||_{\infty}<\infty$ where $||K||_{\infty}=\sup_{x}|K(x)|$. The KDE for $f_{X}$ is $$\begin{aligned}
\ft X(\mathbf{X}_{j}) & = & \frac{1}{Mh_{X}^{d_{X}}}\sum_{\substack{i=1\\
i\neq j
}
}^{N}K_{X}\left(\frac{\mathbf{X}_{j}-\mathbf{X}_{i}}{h_{X}}\right).\label{eq:fx}\end{aligned}$$ The KDEs $\ft Y(\Y_{j})$ and $\ft Z(\mathbf{X}_{j},\mathbf{Y}_{j})$ (where $h_{Z}=(h_{X},h_{Y})$) for estimating $f_{Y}$ and $f_{XY}$, respectively, are defined similarly using $K_{Y}$ and the product kernel $K_{X}\cdot K_{Y}$. Then $G_{1}(\mathbf{X};\mathbf{Y})$ is estimated as $$\gt=\frac{1}{N}\sum_{i=1}^{N}g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right).\label{eq:Gest}$$

Convergence Rates
-----------------

To derive the convergence rates of $\gt$ we assume that 1) $f_{X}$, $f_{Y}$, $f_{XY}$, and $g$ are smooth; 2) $f_{X}$ and $f_{Y}$ have bounded support sets $\mathcal{S}_{X}$ and $\mathcal{S}_{Y}$; 3) $f_{X}$, $f_{Y}$, and $f_{XY}$ are strictly lower bounded on their support sets. More specifically, we assume that the densities belong to the bounded Hölder class $\Sigma(s,H)$ (the precise definition is included in the appendices) which implies that the densities are $r=\left\lfloor s\right\rfloor $ times differentiable. These assumptions are comparable to those in similar studies on asymptotic convergence analysis [@moon2016isit; @moon2014nips; @moon2014isit; @singh2014renyi; @singh2014exponential; @sricharan2013ensemble; @krishnamurthy2014divergence; @kandasamy2015nonparametric]. To derive the convergence rates without boundary corrections, we also assume that 4) the boundary of the support set is smooth with respect to the corresponding kernels. The full assumptions are

-   $(\mathcal{A}.0)$: The kernels $K_{X}$ and $K_{Y}$ are symmetric product kernels with bounded support.

-   $(\mathcal{A}.1)$: There exist constants $\epsilon_{0},\epsilon_{\infty}$ such that $0<\epsilon_{0}\leq f_{X}(x)\leq\epsilon_{\infty}<\infty,\,\forall x\in\mathcal{S}_{X}$, $\epsilon_{0}\leq f_{Y}(y)\leq\epsilon_{\infty},\,\forall y\in\mathcal{S}_{Y}$, and $\epsilon_{0}\leq f_{XY}(x,y)\leq\epsilon_{\infty},\,\forall(x,y)\in\mathcal{S}_{X}\times\mathcal{S}_{Y}$.

-   $(\mathcal{A}.2)$: Each of the densities belong to $\Sigma(s,H)$ in the interior of their support sets with $s\geq2$.

-   $(\mathcal{A}.3)$: $g\left(t_{1}/t_{2}\right)$ has an infinite number of mixed derivatives wrt $t_{1}$ and $t_{2}$.

-   $(\mathcal{A}.4$): $\left|\frac{\partial^{k+l}g(t_{1},t_{2})}{\partial t_{1}^{k}\partial t_{2}^{l}}\right|/(k!l!)$, $k,l=0,1,\ldots$ are strictly upper bounded for $\epsilon_{0}\leq t_{1},t_{2}\leq\epsilon_{\infty}$.

-   $(\mathcal{A}.5)$: Let $K$ be either $K_{X}$ or $K_{Y}$, $\mathcal{S}$ either $\mathcal{S}_{X}$ or $\mathcal{S}_{Y}$, $h$ either $h_{X}$ or $h_{Y}$. Let $p_{x}(u):\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a polynomial in $u$ of order $q\leq r=\left\lfloor s\right\rfloor $ whose coefficients are a function of $x$ and are $r-q$ times differentiable. For any positive integer $t$ $$\int_{x\in\mathcal{S}}\left(\int_{u:K(u)>0,\,x+uh\notin\mathcal{S}}K(u)p_{x}(u)du\right)^{t}dx=v_{t}(h),$$ where $v_{t}(h)$ admits the expansion $$v_{t}(h)=\sum_{i=1}^{r-q}e_{i,q,t}h^{i}+o\left(h^{r-q}\right),$$ for some constants $e_{i,q,t}$.

Assumption $(\mathcal{A}.5)$ states that the support of the density is smooth with respect to the kernel $K$ in the sense that the expectation with respect to any random variable $u$ of the area of the kernel that falls outside the support $\mathcal{S}$ is a smooth function of the bandwidth $h$ provided that the distribution function $p_{x}(u)$ of $u$ is smooth (e.g. $s\geq2$). The inner integral 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 $p_{x}(u)$. As an example, this smoothness assumption is satisfied when the support is rectangular and the kernel is the uniform rectangular kernel [@moon2016arxiv; @moon2016isit]. Note that this boundary assumption does not result in parametric convergence rates for the plug-in estimator $\gt$, which is in contrast with the boundary assumptions in [@singh2014exponential; @singh2014renyi; @krishnamurthy2014divergence; @kandasamy2015nonparametric]. However, the estimators in [@singh2014exponential; @singh2014renyi; @krishnamurthy2014divergence; @kandasamy2015nonparametric] 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 $\gt$ without boundary correction.

\[thm:bias\](Bias) Under assumptions $\mathcal{A}.0-\mathcal{A}.5$ and for general $g$, the bias of $\gt$ is $$\begin{aligned}
\bias\left[\gt\right] & = & \sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{10,i,j}h_{X}^{i}h_{Y}^{j}+\frac{c_{11}}{Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}}\nonumber \\
 &  & +O\left(h_{X}^{s}+h_{Y}^{s}+\frac{1}{Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}}\right).\label{eq:bias1}\end{aligned}$$ If $g\left(t_{1},t_{2}\right)$ also has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$ for each $1\leq k,l,\leq\lambda$, the bias of $\gt$ is $$\begin{aligned}
 &  & \lefteqn{\bias\left[\gt\right]}\nonumber \\
 &  & =\sum_{\substack{m,n=0\\
i+j+m+n\neq0
}
}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i,j=0}^{r}c_{11,j,i,m,n}\frac{h_{X}^{i}h_{Y}^{j}}{\left(Nh_{X}^{d_{X}}\right)^{m}\left(Nh_{Y}^{d_{Y}}\right)^{n}}\nonumber \\
 &  & +\sum_{m=1}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i=0}^{r}\sum_{j=0}^{r}c_{13,m,n,j}h_{X}^{i}h_{Y}^{j}/\left(Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\right)^{m}\nonumber \\
 &  & +O\left(h_{X}^{s}+h_{Y}^{s}+1/\left(Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\right)^{\lambda/2}\right).\label{eq:bias2}\end{aligned}$$

The constants in both (\[eq:bias1\]) and (\[eq:bias2\]) depend only on the densities and their derivatives, the functional $g$ and its derivatives, and the kernels. They are independent of $N,$ $h_{X}$, and $h_{Y}.$

The purpose of Theorem \[thm:bias\] is two-fold. First, we use Theorem \[thm:bias\] to derive the bias expressions for the MI plug-in estimators when $\mathbf{X}$ and $\mathbf{Y}$ may have a mixture of discrete and continuous components (cases 2 and 3) in Section \[sec:mixed\]. Second, in conjunction with Theorem \[thm:variance\] which follows, the results in Theorem \[thm:bias\] can be used to derive MI ensemble estimators in Appendix \[subsec:cont\_ensemble\] that achieve the parametric MSE convergence rate when the densities are sufficiently smooth. The expression in (\[eq:bias2\]) enables us to achieve the parametric rate under less restrictive smoothness assumptions on the densities ($s>d/2$ for (\[eq:bias2\]) compared to $s\geq d$ for (\[eq:bias1\])). The extra condition required on the mixed derivatives of $g$ to obtain the expression in (\[eq:bias2\]) is satisfied, for example, for Shannon and Renyi information measures.

\[thm:variance\](Variance) If the functional $g$ is Lipschitz continuous in both of its arguments with Lipschitz constant $C_{g}$, then the variance of $\gt$ is $$\var\left[\gt\right]\leq\frac{22C_{g}^{2}||K_{X}\cdot K_{Y}||_{\infty}^{2}}{N}.$$

Similar to Theorem \[thm:bias\], Theorem \[thm:variance\] is used to derive variance expressions for the MI plug-in estimators under cases 2 and 3. Theorem \[thm:variance\] is also necessary to derive optimally weighted ensemble estimators. The proofs of Theorems \[thm:bias\] and \[thm:variance\] are similar to the proofs of the bias and variance results for the divergence functional estimators in [@moon2016arxiv]. The primary difference is in handling certain products of the marginal KDEs that appear in the expansion of the MSE. See Appendix \[sec:biasProof\] and \[sec:VarProof\] for details.

Theorems \[thm:bias\] and \[thm:variance\] indicate that for the MSE of the plug-in estimator to go to zero for case 1, we require $h_{X},h_{Y}\rightarrow0$ and $Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\rightarrow\infty$. The Lipschitz assumption on $g$ is comparable to other nonparametric estimators of distributional functionals [@kandasamy2015nonparametric; @singh2014exponential; @singh2014renyi; @moon2016arxiv; @krishnamurthy2014divergence]. Specifically, assumption $\mathcal{A}.1$ ensures that functionals such as those for Shannon and Renyi informations are Lipschitz on the space $\epsilon_{0}$ to $\epsilon_{\infty}$.

Mixed Random Variables {#sec:mixed}
======================

In this section, we extend the results of Section \[sec:MI\_est\] to MI estimation when $\mathbf{X}$ and $\mathbf{Y}$ may have a mixture of discrete and continuous components. For simplicity, we focus primarily on the important case when $\mathbf{X}$ is continuous and $\mathbf{Y}$ is discrete (case 2 in Section \[sec:intro\]). The more general case when $\X$ and $\Y$ may have any mixture of continuous and discrete components (case 3 in Section \[sec:intro\]) is discussed in Section \[subsec:general\_rates\]. As an example of the former case, if $\mathbf{Y}$ is a predictor variable (e.g. classification labels), then the MI between $\mathbf{X}$ and $\mathbf{Y}$ indicates the value of $\mathbf{X}$ as a predictor of $\mathbf{Y}$. Although $\mathbf{Y}$ is discrete, $f_{XY}=f_{Z}$ is also a density. Let $\mathcal{S}_{X}$ be the support of the density $f_{X}$ and $\mathcal{S}_{Y}$ be the support of the probability mass function $f_{Y}$. The MI is $$\begin{aligned}
 & \lefteqn{G_{2}\left(\mathbf{X};\mathbf{Y}\right)}\nonumber \\
 & = & \sum_{y\in\mathcal{S}_{Y}}\int g\left(\frac{f_{X}(x)f_{Y}(y)}{f_{XY}(x,y)}\right)f_{XY}(x,y)dx\label{eq:MI_cond}\\
 & = & \sum_{y\in\mathcal{S}_{Y}}f_{Y}(y)\int g\left(\frac{f_{X}(x)}{f_{X|Y}(x|y)}\right)f_{X|Y}(x|y)dx.\nonumber \end{aligned}$$

Let $\mathbf{N}_{y}=\sum_{i=1}^{N}1_{\left\{ \mathbf{Y}_{i}=y\right\} }$ where $y\in\mathcal{S}_{Y}$. Let $\ft X$ be as in (\[eq:fx\]) and define $\mathcal{X}_{y}=\left\{ \mathbf{X}_{i}\in\left\{ \mathbf{X}_{1},\dots,\mathbf{X}_{N}\right\} |\mathbf{Y}_{i}=y\right\} $. Then if $\mathbf{X}_{i}\in\mathcal{X}_{y}$, the KDE of $f_{X|Y}(x|y)$ is $$\begin{aligned}
\ft{X|y}(\mathbf{X}_{i}) & = & \frac{1}{\left(\mathbf{N}_{y}-1\right)h_{X|y}^{d_{X}}}\sum_{\substack{\mathbf{X}_{j}\in\mathcal{X}_{y}\\
i\neq j
}
}K_{X}\left(\frac{\mathbf{X}_{i}-\mathbf{X}_{j}}{h_{X|y}}\right).\end{aligned}$$ We define the plug-in estimator $\g{h_{X},h_{X|Y}}$ of (\[eq:MI\_cond\]) as $$\begin{aligned}
\g{h_{X},h_{X|y}} & =\frac{1}{\mathbf{N}_{y}}\sum_{\mathbf{X}\in\mathcal{X}_{y}}g\left(\ft X(\mathbf{X})/\ft{X|y}(\mathbf{X})\right),\nonumber \\
\implies\g{h_{X},h_{X|Y}} & =\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\g{h_{X},h_{X|y}}.\label{eq:mixed_est}\end{aligned}$$

Convergence Rates {#subsec:mixed_conv}
-----------------

To apply the theory of optimally weighted ensemble estimation to $\g{h_{X},h_{X|Y}}$, we need to know its MSE as a function of the bandwidths and the sample size.

\[thm:bias\_mixed\](Bias) Assume that assumptions $\mathcal{A}.0-\mathcal{A}.5$ apply to the functional $g$, the kernel $K_{X}$, and the densities $f_{X}$ and $f_{X|Y}$. Assume that $\mathbf{h}_{X|y}=l\mathbf{N}_{y}^{-\beta}$ with $0<\beta<\frac{1}{d_{X}}$ and $l$ a positive number. Then the bias of $\g{h_{X},h_{X|Y}}$ is $$\begin{aligned}
 & \lefteqn{\bias\left[\g{h_{X},h_{X|Y}}\right]}\nonumber \\
 & =\sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{13,i,j}h_{X}^{i}l^{j}N^{-j\beta}+\frac{c_{14,X}}{Nh_{X}^{d_{X}}}+\frac{c_{14,Y}}{l^{d_{X}}N^{1-\beta d_{X}}}\nonumber \\
 & +O\left(h_{X}^{s}+N^{-s\beta}+\frac{1}{Nh_{X}^{d_{X}}}+\frac{1}{N^{1-\beta d_{X}}}\right).\label{eq:bias_mixed1}\end{aligned}$$ If $g\left(t_{1},t_{2}\right)$ also has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$ for each $1\leq j,l\leq\lambda$, then the bias is $$\begin{aligned}
 & \lefteqn{\bias\left[\g{h_{X},h_{X|Y}}\right]}\nonumber \\
 & =\sum_{\substack{m,n=0\\
i+j+m+n\neq0
}
}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i,j=0}^{r}c_{14,j,i,m,n}\frac{h_{X}^{i}l^{j}N^{-j\beta}}{\left(Nh_{X}^{d_{X}}\right)^{m}\left(l^{d_{X}}N^{1-\beta d_{X}}\right)^{n}}\nonumber \\
 & +O\left(h_{X}^{s}+N^{-s\beta}+\frac{1}{\left(Nh_{X}^{d_{X}}\right)^{\lambda/2}}+\frac{1}{\left(N^{1-\beta d_{X}}\right)^{\lambda/2}}\right).\label{eq:bias_mixed2}\end{aligned}$$

We focus on (\[eq:bias\_mixed1\]) as (\[eq:bias\_mixed2\]) follows similarly. It can be shown that $$\bias\left[\g{h_{X},h_{X|Y}}\right]=\bE\left[\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\bias\left[\left.\g{h_{X},h_{X|y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right].$$ The conditional bias of $\g{h_{X},h_{X|y}}$ given $\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}$ can then be obtained from Theorem \[thm:bias\] as $$\begin{aligned}
 & \lefteqn{\bias\left[\left.\g{h_{X},h_{X|y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]}\\
 & =\sum_{\substack{i,j=0\\
i+j\neq0
}
}^{r}c_{10,i,j}h_{X}^{i}\mathbf{h}_{X|y}^{j}\\
 & +O\left(h_{X}^{s}+\mathbf{h}_{X|y}^{s}+\frac{1}{\mathbf{N}_{y}h_{X}^{d_{X}}}+\frac{1}{\mathbf{N}_{y}\mathbf{h}_{X|y}^{d_{X}}}\right)\end{aligned}$$ Then given that $\mathbf{h}_{X|y}\propto\mathbf{N}_{y}^{-\beta}$, (\[eq:mixed\_est\]) gives terms of the form of $\mathbf{N}_{y}^{1-\gamma}$ with $\gamma>0$. $\mathbf{N}_{y}$ is a binomial random variable with parameter $f_{Y}(y)$, $N$ trials, and mean $Nf_{Y}(y)$. Thus we need to compute the fractional moments of a binomial random variable. By the generalized binomial theorem, we have that $$\begin{aligned}
\mathbf{N}_{y}^{\alpha} & =\left(\mathbf{N}_{y}-Nf_{Y}(y)+Nf_{Y}(y)\right)^{\alpha}\nonumber \\
 & =\sum_{i=0}^{\infty}\left(\begin{array}{c}
\alpha\\
i
\end{array}\right)\left(Nf_{Y}(y)\right)^{\alpha-i}\left(\mathbf{N}_{y}-Nf_{Y}(y)\right)^{i},\nonumber \\
 & \lefteqn{\implies\bE\left[\mathbf{N}_{y}^{\alpha}\right]}\nonumber \\
 & =\sum_{i=0}^{\infty}\left(\begin{array}{c}
\alpha\\
i
\end{array}\right)\left(Nf_{Y}(y)\right)^{\alpha-i}\bE\left[\left(\mathbf{N}_{y}-Nf_{Y}(y)\right)^{i}\right].\label{eq:fractional_moment}\end{aligned}$$ From [@riordan1937moment], the $i$-th central moment of $\mathbf{N}_{y}$ has the form of $$\bE\left[\left(\mathbf{N}_{Y}-Nf_{Y}(y)\right)^{i}\right]=\sum_{n=0}^{\left\lfloor i/2\right\rfloor }c_{n,i}(f_{Y}(y))N^{n}.$$ Thus $\bE\left[\mathbf{N}_{y}^{1-\gamma}\right]$ has terms proportional to $N^{1-\gamma-i+n}\leq N^{1-\gamma-\left\lfloor i/2\right\rfloor }$ for $i=0,1,\dots$ since $n\leq\left\lfloor i/2\right\rfloor $. Then since there is an $N$ in the denominator of (\[eq:mixed\_est\]), this leaves terms of the form of $N^{-\gamma}$ when $i=0,1$ and $N^{-1}$ for $i\geq2$. This completes the proof for the bias. See Appendix \[sec:MixedProofs\] for more details.

\[thm:var\_mixed\]If the functional $g$ is Lipschitz continuous in both of its arguments and $\mathcal{S}_{Y}$ is finite, then the variance of $\g{h_{X},h_{X|Y}}$ is $O(1/N)$.

By the law of total variance, we have $$\begin{aligned}
\var\left[\g{h_{X},h_{X|Y}}\right] & =\bE\left[\var\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right]\\
 & +\var\left[\bE\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right].\end{aligned}$$ Given all of the $\mathbf{Y}_{i}$’s, the estimators $\g{h_{X},h_{X|y}}$ are all independent since they use different sets of $\mathbf{X}_{i}$’s for each $y$. From Theorem \[thm:variance\], we know that $\var\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]=O\left(\sum_{y\in\mathcal{S}_{Y}}\mathbf{N}_{y}/N^{2}\right)$. Taking the expectation then yields $O(1/N)$.

For the second term, we know from the proof of Theorem \[thm:bias\_mixed\] that $\bE\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]$ yields a sum of terms of the form of $\mathbf{N}_{y}^{\gamma}/N$ for $0<\gamma\leq1$. Taking the variance of the sum of these terms yields a sum of terms of the form $\var\left[\mathbf{N}_{y}^{\gamma}\right]/N^{2}$ (the covariance terms can be bounded by the Cauchy-Schwarz inequality to yield similar terms). Then $\var\left[\mathbf{N}_{y}^{\gamma}\right]$ can be bounded by taking a Taylor series expansion of the functions $\mathbf{N}_{y}^{\gamma}$ and $\mathbf{N}_{y}^{2\gamma}$ at the point $Nf_{Y}(y)$ which yields an expression that depends on the central moments of $\mathbf{N}_{y}$. From this, we obtain $\var\left[\mathbf{N}_{y}^{\gamma}\right]=O(N)$ which completes the proof. See Appendix \[sec:MixedProofs\] for details.

Theorems \[thm:bias\_mixed\] and \[thm:var\_mixed\] provide exact expressions for the bias and bounds on the variance of the plug-in MI estimator, respectively. It is shown in Section \[sec:mixed\_ensemble\] that the MSE of the plug-in estimator converges very slowly to zero under this setting. However, Theorems \[thm:bias\_mixed\] and \[thm:var\_mixed\] provide with us the necessary information for applying the theory of optimally weighted ensemble estimation to obtain estimators with improved rates. This is done in Section \[sec:mixed\_ensemble\].

Extension to Other Cases {#subsec:general_rates}
------------------------

The results in Section \[subsec:mixed\_conv\] can be extended to the case where $\mathbf{X}$ and/or $\mathbf{Y}$ may have a mixture of continuous and discrete components (case 3 in Section \[sec:intro\]). This scenario can be divided further into three different cases: A) $\X$ is continuous and $\Y$ has a mixture of discrete and continuous components; B) $\X$ and $\Y$ both have a mixture of discrete and continuous components; C) $\Y$ is discrete and $\X$ has a mixture of discrete and continuous components. Consider case A first. Denote the discrete and continuous components of $\mathbf{Y}$ as $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$, respectively. Denote the respective support sets as $\mathcal{S}_{Y_{1}}$ and $\mathcal{S}_{Y_{2}}$. We can then write $$\begin{aligned}
 & \lefteqn{G_{3A}(\mathbf{X};\mathbf{Y})}\nonumber \\
 & =\sum_{y_{1}\in\mathcal{S}_{Y_{1}}}\int g\left(\frac{f_{X}(x)f_{Y}(y_{1},y_{2})}{f_{XY}(x,y_{1},y_{2})}\right)f_{XY}(x,y_{1},y_{2})dxdy_{2}\nonumber \\
 & =\sum_{y_{1}\in\mathcal{S}_{Y_{1}}}f_{Y_{1}}(y_{1})\int g\left(\frac{f_{X}(x)f_{Y_{2}|Y_{1}}(y_{2}|y_{1})}{f_{XY_{2}|Y_{1}}(x,y_{2}|y_{1})}\right)\nonumber \\
 & \times f_{XY_{2}|Y_{1}}(x,y_{2}|y_{1})dxdy_{2}.\label{eq:mixed_general}\end{aligned}$$ The subscript $3A$ indicates that we are considering case A under the third case described in the introduction. The expression in (\[eq:mixed\_general\]) is very similar to the expression in (\[eq:MI\_cond\]). After plugging in KDEs for the corresponding densities and conditional densities, a nearly identical procedure to that in Section \[subsec:mixed\_conv\] can be followed to derive the bias and variance of the corresponding plug-in estimator.

Now consider case B. Denote the discrete and continuous components of $\mathbf{X}$ as $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, respectively. Then if $\mathbf{Y}_{1}$ is the discrete component of $\mathbf{Y}$, then the expression inside the $g$ functional in (\[eq:mixed\_general\]) includes $f_{X_{1}}(x_{1})f_{Y_{1}}(y_{1})/f_{X_{1}Y_{1}}(x_{1},y_{1})$. Thus the plug-in estimator must include estimators for $f_{X_{1}}(x_{1}),$ $f_{Y_{1}}(y_{1})$, and $f_{X_{1}Y_{1}}(x_{1},y_{1})$. Define $\mathbf{N}_{y_{1}}=\sum_{i=1}^{N}1_{\{\mathbf{Y}_{1,i}=y_{1}\}}$ where $\mathbf{Y}_{1,i}$ is the discrete component of $\mathbf{Y}_{i}$. Then the estimator we use for $f_{Y_{1}}(y_{1})$ is $\mathbf{N}_{y_{1}}/N$. The estimators for $f_{X_{1}}(x_{1})$ and $f_{X_{1}Y_{1}}(x_{1},y_{1})$ are defined similarly. The bias and variance expressions of this plug-in estimator can then be derived with some slight modifications of Theorems \[thm:bias\] and \[thm:variance\]. See Appendix \[subsec:generalCase\] for an expression for $G_{3B}(\mathbf{X};\mathbf{Y})$ in this case and a sketch of these modifications. Case C follows similarly as the expression inside the $g$ functional in (\[eq:mixed\_general\]) includes $f_{X_{1}}(x_{1})f_{Y}(y)/f_{X_{1}Y}(x_{1},y)$ where all the terms are probability mass functions.

The resulting bias and variance expressions in these settings are analogous to those in Theorems \[thm:bias\], \[thm:variance\], and \[thm:bias\_mixed\] as the variance will be $O(1/N)$ and the bias will depend on expansions of the bandwidths for the various KDEs. Ensemble methods can then be applied to improve the MSE convergence rates as described in the next section.

Ensemble Estimation of MI\[sec:mixed\_ensemble\]
================================================

Mixed Random Variables {#subsec:mixed_ensemble}
----------------------

We again focus on the case where $\mathbf{X}$ is continuous and $\mathbf{Y}$ is discrete (case 2 in Section \[sec:intro\]). If no bias correction is performed, then Theorem \[thm:bias\_mixed\] shows that the optimal bias rate of the plug-in estimator $\g{h_{X},h_{X|Y}}$ is $O\left(1/N^{1/(d_{X}+1)}\right)$, which converges very slowly to zero when $d_{X}$ is not small. We use the theory of optimally weighted ensemble estimation to improve this rate. An ensemble of estimators is formed by choosing different bandwidth values. Consider first the case where (\[eq:bias\_mixed1\]) applies. Let $\mathcal{L}$ be a set of real positive numbers with $|\mathcal{L}|=L$. This set will parameterize the bandwidths for $\ft X$ and $\ft{X|y}$ resulting in $L$ estimators in the ensemble. While different parameter sets for $\ft X$ and $\ft{X|y}$ can be chosen, we only use one set here for simplicity of exposition. To ensure that the final terms in (\[eq:bias\_mixed1\]) are $O(1/\sqrt{N})$ when $s\geq d$, for each estimator in the ensemble we choose $h_{X}(l)=lN^{-1/(2d_{X})}$ and $\mathbf{h}_{X|y}(l)=l\mathbf{N}_{y}^{-1/(2d_{X})}$ where $l\in\mathcal{L}$. Define $w$ to be a weight vector parameterized by $l\in\mathcal{L}$ with $\sum_{l\in\mathcal{L}}w(l)=1$ and define $$\g{w,1}=\sum_{l\in\mathcal{L}}w(l)\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\g{h_{X}(l),h_{X|y}(l)}.\label{eq:ensemble}$$ From Theorem \[thm:bias\_mixed\], the bias of $\g{w,1}$ is $$\begin{aligned}
\bias\left[\g{w,1}\right] & =\sum_{l\in\mathcal{L}}\sum_{i=1}^{r}\theta\left(w(l)l^{i}N^{\frac{-i}{2d_{X}}}\right)\nonumber \\
 & +O\left(\sqrt{L}||w||_{2}\left(N^{\frac{-s}{2d_{X}}}+N^{\frac{-1}{2}}\right)\right),\label{eq:weight_bias}\end{aligned}$$ where we use $\theta$ notation to omit the constants.

We use the general theory of optimally weighted ensemble estimation in [@moon2016isit] to improve the MSE convergence rate of the plug-in estimator by using the weights to cancel the lower order terms in (\[eq:weight\_bias\]). The theory is as follows. Let $\left\{ \hat{\mathbf{E}}_{l}\right\} _{l\in\mathcal{L}}$ be an indexed ensemble of estimators with the weighted ensemble estimator $\hat{\mathbf{E}}_{w}=\sum_{l\in\mathcal{L}}w(l)\hat{\mathbf{E}}_{l}$ satisfying:

-   $\mathcal{C}.1$. Let $c_{i}$ be constants depending on the underlying density, $J=\{i_{1},\dots i_{I}\}$ a finite index set with $I<L$, $\psi_{i}(l)$ basis functions depending only on the parameter $l$ and not on $N$, $\phi_{i}(N)$ functions of the sample size $N$ that are independent of $l$. Assume the bias is $$\bias\left[\hat{\mathbf{E}}_{l}\right]=\sum_{i\in J}c_{i}\psi_{i}(l)\phi_{i}(N)+O\left(\frac{1}{\sqrt{N}}\right).$$

-   $\mathcal{C}.2$. Assume the variance is $$\var\left[\hat{\mathbf{E}}_{l}\right]=c_{v}\left(\frac{1}{N}\right)+o\left(\frac{1}{N}\right).$$

[@moon2016isit] \[thm:opt\_weight\]If conditions $\mathcal{C}.1$ and $\mathcal{C}.2$ hold for an ensemble of estimators $\left\{ \hat{\mathbf{E}}_{l}\right\} _{l\in\mathcal{L}}$, then there exists a weight vector $w_{0}$ such that the MSE of $\hat{\mathbf{E}}_{w_{0}}$ attains the parametric rate of convergence of $O\left(1/N\right)$. The weight $w_{0}$ is the solution to the offline convex optimization problem *$$\begin{array}{rl}
\min_{w} & ||w||_{2}\\
subject\,to & \sum_{l\in\mathcal{L}}w(l)=1,\\
 & \gamma_{w}(i)=\sum_{l\in\mathcal{L}}w(l)\psi_{i}(l)=0,\,i\in J.
\end{array}\label{eq:optimize}$$*

To apply Theorem \[thm:opt\_weight\] to an ensemble of estimators, all $\phi_{i}(N)$ functions that converge to zero slower than $1/\sqrt{N}$ and the corresponding $\psi_{i}(l)$ functions must be known for the base estimator. Otherwise, Theorem \[thm:opt\_weight\] can only be guaranteed to improve the bias up to the slowest unknown bias rate. This theorem was applied in [@moon2016isit] to the problem of divergence functional estimation where the plug-in estimator has slowly converging bias but the resulting ensemble estimator achieves the parametric rate for sufficiently smooth densities.

We apply Theorem \[thm:opt\_weight\] to the ensemble estimator $\g{w,1}$ as conditions $\mathcal{C}.1$ and $\mathcal{C}.2$ are satisfied with $\phi_{i}(N)=N^{-i/(2d_{X})}$ and $\psi_{i}(l)=l^{i}$ for $i\in\{1,\dots r\}$ as seen in (\[eq:bias\_mixed1\]) and (\[eq:weight\_bias\]). If $s\geq d_{X}$, then the MSE of the optimally weighted estimator $\g{w_{0},1}$ is $O(1/N)$. A similar approach can be used for the case where $\X$ contains a mixture of continuous and discrete components and $\Y$ is discrete (or vice versa). To the best of our knowledge, these are the first nonparametric estimators to achieve the MSE parametric rate in this setting of mixed random variables.

If the mixed derivatives of the functional $g$ satisfy the extra condition required for (\[eq:bias\_mixed2\]), we can define an ensemble estimator $\g{w_{0},2}$ that achieves the parametric MSE rate if $s>d_{X}/2$. For simplicity, we focus primarily on $\g{w_{0},1}$. See Appendix \[subsec:Odin2\] for details on $\g{w_{0},2}$.

In practice, the optimization problem in (\[eq:optimize\]) typically results in a very large increase in variance. Thus we follow the lead of [@moon2016arxiv; @moon2014isit; @moon2014nips; @sricharan2013ensemble] and use a relaxed version of (\[eq:optimize\]): $$\begin{array}{rl}
\min_{w} & \epsilon\\
subject\,to & \sum_{l\in\mathcal{L}}w(l)=1,\\
 & \left|\gamma_{w}(i)N^{\frac{1}{2}}\phi_{i}(N)\right|\leq\epsilon,\,\,i\in J,\\
 & \left\Vert w\right\Vert _{2}^{2}\leq\eta.
\end{array}\label{eq:relaxed}$$ As shown in [@moon2016arxiv; @moon2014isit; @moon2014nips; @sricharan2013ensemble], the ensemble estimator $\g{w_{0},1}$ using the resulting weight vector from the optimization problem in (\[eq:relaxed\]) still achieves the parametric MSE convergence rate under the same assumptions as described previously. It was also shown in [@moon2016arxiv] that the heuristic of setting $\eta=\epsilon$ works well in practice. Algorithm \[alg:estimator\] summarizes the estimator $\g{w_{0},1}$.

$L$ positive real numbers $\mathcal{L}$, samples $\left\{ \mathbf{Z}_{1},\dots,\mathbf{Z}_{N}\right\} $ from $f_{XY}$, dimension $d_{X}$, function $g$, kernel $K_{X}$

The optimally weighted MI estimator $\g{w_{0},1}$

Solve for $w_{0}$ using (\[eq:relaxed\]) with basis functions $\psi_{i}(l)=l^{i}$, $\phi_{i}(N)=N^{-i/(2d_{X})},$ $l\in\mathcal{L}$, and $0\leq i\leq d_{X}$.

$\mathbf{N}_{y}\leftarrow\sum_{i=1}^{N}1_{\{\mathbf{Y}_{i}=y\}}$

$h_{X}(l)\leftarrow lN^{-1/(2d_{X})},$ $\mathbf{h}_{X|y}(l)\leftarrow l\mathbf{N}_{y}^{-1/(2d_{Y})}$

Calculate $\ftl X(\mathbf{X}_{i})$, $\ftl{X|y}(\mathbf{X}_{i})$ as described in the text

$\g{h_{X}(l),\mathbf{h}_{X|y}(l)}\leftarrow\frac{1}{\mathbf{N}_{y}}\sum_{\mathbf{X}\in\mathcal{X}_{y}}g\left(\frac{\ftl X(\mathbf{X})}{\ftl{X|y}(\mathbf{X})}\right)$

$\g{w_{0},1}\leftarrow\sum_{l\in\mathcal{L}}w_{0}(l)\sum_{y\in\mathcal{S}_{y}}\frac{\mathbf{N}_{y}}{N}\g{h_{X}(l),\mathbf{h}_{X|y}(l)}$

A similar approach can be used to derive an ensemble estimator $\g{w_{0},1}^{cont}$ for the case when $\X$ and $\Y$ are continuous (case 1 in Section \[sec:intro\]). See Appendix \[subsec:cont\_ensemble\] for details. The case where $\X$ and $\Y$ both contain a mixture of discrete and continuous components follows similarly.

Parameter Selection
-------------------

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

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

2.  Ensure the bandwidths are sufficiently distinct. Similar bandwidth values lead to negligible decrease in bias and many bandwidth values may increase $||w_{0}||_{2}$ resulting in an increase in variance [@sricharan2013ensemble].

3.  Select $L=|\mathcal{L}|>|J|=I$ to obtain a feasible solution for the optimization problems in (\[eq:optimize\]) and (\[eq:relaxed\]). We find that choosing a value of $30\leq L\leq60$, and setting $\mathcal{L}$ to be $L$ linearly spaced values between the minimum and maximum values described above works well in practice.

The resulting ensemble estimators are robust in the sense that they are not sensitive to the exact choice of the bandwidths or the number of estimators as long as the the rough rules-of-thumb given above are followed. Moon et al [@moon2016arxiv; @moon2016isit] gives more details on ensemble estimator parameter selection for continuous divergence estimation. These details also apply to the continuous parts of the mixed cases for MI estimation in this paper.

Since the optimal weight $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\left(N^{2}\right)$ using the standard implementation. For very large datasets, more efficient KDE implementations (e.g. [@raykar2010fast]) can be used to reduce the computational burden.

Central Limit Theorem {#subsec:clt}
---------------------

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

\[thm:clt\] Let $\g w^{cont}$ be a weighted ensemble estimator when $\X$ and $\Y$ are continuous with bandwidths $h_{X}(l_{X})$ and $h_{Y}(l_{Y})$ for each estimator in the ensemble. Assume that the functional $g$ is Lipschitz in both arguments with Lipschitz constant $C_{g}$ and that $h_{X}(l_{X}),\,h_{Y}(l_{Y})=o(1)$, $N\rightarrow\infty$, and $Nh_{X}^{d_{X}}(l_{X}),\,Nh_{Y}^{d_{Y}}(l_{Y})\rightarrow\infty$ for each $l_{X}\in\mathcal{L}_{X}$ and $l_{Y}\in\mathcal{L}_{Y}$. Then for fixed $\mathcal{L}_{X}$ and $\mathcal{L}_{Y}$, and if $\mathbf{S}$ is a standard normal random variable, $$\Pr\left(\left(\g w^{cont}-\bE\left[\g w^{cont}\right]\right)/\sqrt{\var\left[\g w^{cont}\right]}\leq t\right)\rightarrow\Pr\left(\mathbf{S}\leq t\right).$$

The proof is based on an application of Slutsky’s Theorem preceded by an application of the Efron-Stein inequality (see Appendix \[sec:cltProof\]).

If the space $\mathcal{S}_{Y}$ is finite, then the ensemble estimators for the mixed component case also obey a central limit theorem. The proof follows by an application of Slutsky’s Theorem combined with Theorem \[thm:clt\].

Let $\g w$ be a weighted ensemble estimator when $\X$ is continuous and $\Y$ is discrete with bandwidths $h_{X}(l)$ and $h_{X|y}(l)$ for each estimator in the ensemble. Assume that the functional $g$ is Lipschitz in both arguments and that $h_{X},\,h_{X|y}=o(1)$, $N\rightarrow\infty$, and $Nh_{X}^{d_{X}},\,Nh_{X|y}^{d_{X}}\rightarrow\infty$ for each $l\in\mathcal{L}$ and $\forall y\in\mathcal{S}_{Y}$ with $\mathcal{S}_{y}$ finite. Then for fixed $\mathcal{L}$, $$\Pr\left(\left(\g w-\bE\left[\g w\right]\right)/\sqrt{\var\left[\g w\right]}\leq t\right)\rightarrow\Pr\left(\mathbf{S}\leq t\right).$$

Experimental Validation {#sec:experiments}
=======================

In this section, we validate our theory by estimating the Rényi-$\alpha$ MI integral (i.e. $g(x)=x^{\alpha}$ in (\[eq:MI\_cond\]); see [@principe2010information]) where $\mathbf{X}$ is a mixture of truncated Gaussian random variables restricted to the unit cube and $\mathbf{Y}$ is a categorical random variable. We choose Rényi MI as it has received recent interest (e.g. [@pal2010estimation]) and the estimation problem does not reduce to entropy estimation in contrast with Shannon MI. Thus this is a clear case where there are no other nonparametric estimators that are known to achieve the parametric MSE rate.

We consider two cases. In the first case, $\mathbf{Y}$ has three possible outcomes (i.e. $|\mathcal{S}_{Y}|=3$) and respective probabilities $\Pr(\mathbf{Y}=0)=\Pr(\mathbf{Y}=1)=2/5$ and $\Pr(\mathbf{Y}=2)=1/5$. The conditional covariance matrices are all $0.1\times I_{d}$ and the conditional means are, respectively, $\bar{\mu}_{0}=0.25\times\bar{1}_{d}$, $\bar{\mu}_{1}=0.75\times\bar{1}_{d}$, and $\bar{\mu}_{2}=0.5\times\bar{1}_{d}$, where $I_{d}$ is the $d\times d$ identity matrix and $\bar{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 $\alpha=0.5$. We set $\mathcal{L}$ to be 40 linearly spaced values between 1.2 and 3. The bandwidth in the KDE plug-in estimator is also set to $2.1N^{-1/(2d)}$.

The top three plots in Figure \[fig:mseplot\] 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 $\g{w_{0},1}$ for various sample sizes and for $d=4,\,6,\,9$, respectively. The ensemble estimator outperforms the standard plug-in estimator, especially for larger sample sizes and larger dimensions. This demonstrates that while an individual kernel estimator performs poorly, an ensemble of estimators including the individual estimator performs well.

For the second case, $\mathbf{Y}$ has six possible outcomes (i.e. $|\mathcal{S}_{Y}|=6$) and respective probabilities $\Pr(\Y=0)=0.35$, $\Pr(\Y=1)=0.2$, $\Pr(\Y=2)=\Pr(\Y=3)=0.15$, $\Pr(\Y=4)=0.1$, and $\Pr(\Y=5)=0.05$. We chose $\alpha=0.5$ and $d=6$. The conditional covariances matrices are again $0.1\times I_{d}$ and the conditional means are, respectively, $\bar{\mu}_{0}=0.25\times\bar{1}_{d}$, $\bar{\mu}_{1}=0.75\times\bar{1}_{d}$, and $\bar{\mu}_{2}=0.5\times\bar{1}_{d}$, $\bar{\mu}_{3}=\left(0.25\times\bar{1}_{4}^{T},0.5\times\bar{1}_{2}^{T}\right)^{T}$, $\bar{\mu}_{4}=\left(0.75\times\bar{1}_{2}^{T},0.375\times\bar{1}_{4}^{T}\right)^{T}$, and $\bar{\mu}_{5}=\left(0.5\times\bar{1}_{4}^{T},0.25\times\bar{1}_{2}^{T}\right)^{T}$. The parameters for the ensemble estimators and the KDE plug-in estimators are the same as in the top three plots in Figure \[fig:mseplot\]. The bottom plot in Figure \[fig:mseplot\] again compares the ensemble estimator to the plug-in KDE estimator. The ensemble estimator also outperforms the plug-in estimator in this setting.

![MSE log-log plots as a function of sample size for the uniform kernel plug-in MI estimator (Kernel) and the proposed optimally weighted ensemble estimator $\protect\g{w_{0},1}$ (Weighted) for the distributions described in the text. The top three plots each correspond to the first case where $|\mathcal{S}_{Y}|=3$ and the bottom plot corresponds to the second case where $|\mathcal{S}_{Y}|=6$. The ensemble estimator outperforms the kernel plug-in estimator, especially for larger sample sizes. Note also that as the dimension increases, the performance gap between the two estimators increases.\[fig:mseplot\]](d4_y3){width="50.00000%"}

![MSE log-log plots as a function of sample size for the uniform kernel plug-in MI estimator (Kernel) and the proposed optimally weighted ensemble estimator $\protect\g{w_{0},1}$ (Weighted) for the distributions described in the text. The top three plots each correspond to the first case where $|\mathcal{S}_{Y}|=3$ and the bottom plot corresponds to the second case where $|\mathcal{S}_{Y}|=6$. The ensemble estimator outperforms the kernel plug-in estimator, especially for larger sample sizes. Note also that as the dimension increases, the performance gap between the two estimators increases.\[fig:mseplot\]](d6_y3){width="50.00000%"}

![MSE log-log plots as a function of sample size for the uniform kernel plug-in MI estimator (Kernel) and the proposed optimally weighted ensemble estimator $\protect\g{w_{0},1}$ (Weighted) for the distributions described in the text. The top three plots each correspond to the first case where $|\mathcal{S}_{Y}|=3$ and the bottom plot corresponds to the second case where $|\mathcal{S}_{Y}|=6$. The ensemble estimator outperforms the kernel plug-in estimator, especially for larger sample sizes. Note also that as the dimension increases, the performance gap between the two estimators increases.\[fig:mseplot\]](d9_y3){width="50.00000%"}

![MSE log-log plots as a function of sample size for the uniform kernel plug-in MI estimator (Kernel) and the proposed optimally weighted ensemble estimator $\protect\g{w_{0},1}$ (Weighted) for the distributions described in the text. The top three plots each correspond to the first case where $|\mathcal{S}_{Y}|=3$ and the bottom plot corresponds to the second case where $|\mathcal{S}_{Y}|=6$. The ensemble estimator outperforms the kernel plug-in estimator, especially for larger sample sizes. Note also that as the dimension increases, the performance gap between the two estimators increases.\[fig:mseplot\]](d6_y6){width="50.00000%"}

Conclusion
==========

We derived the MSE convergence rates for plug-in KDE-based estimators of MI measures between $\mathbf{X}$ and $\mathbf{Y}$ when they have only continuous components and for the case where $\mathbf{Y}$ is discrete and $\mathbf{X}$ is continuous. We also showed how convergence rates can be obtained for the case when $\X$ and/or $\Y$ contain a mixture of discrete and continuous components. Using these rates, we defined ensemble estimators that achieve an MSE rate of $O(1/N)$ when the densities are sufficiently smooth and showed that a central limit theorem also holds. To the best of our knowledge, this is the first nonparametric MI estimator that achieves the MSE convergence rate of $O(1/N)$ in this setting of mixed random variables (i.e. $\X$ and $\Y$ are not both purely discrete or purely continuous).

Hölder Class
============

We derive MSE convergence rates for the plug-in estimators in terms of the smoothness of the densities which we characterize by the Hölder Class.

\[def:holder\]*Let $\mathcal{X}\subset\mathbb{R}^{d}$ be a compact space. For $r=(r_{1},\dots,r_{d}),$ $r_{i}\in\mathbb{N},$ define $|r|=\sum_{i=1}^{d}r_{i}$ and $D^{r}=\frac{\partial^{|r|}}{\partial x_{1}^{r_{1}}\dots\partial x_{d}^{r_{d}}}$. The Hölder class $\Sigma(s,H)$ of functions on $L_{2}(\mathcal{X})$ consists of the functions $f$ that satisfy $$\left|D^{r}f(x)-D^{r}f(y)\right|\leq H\left\Vert x-y\right\Vert ^{s-r},$$ for all $x,\,y\in\mathcal{X}$ and for all $r$ s.t. $|r|\leq\left\lfloor s\right\rfloor $.*

For notation, let $\ez Z$ denote the conditional expectation given $\mathbf{Z}$.

MI Ensemble Estimation Extensions
=================================

Continuous Random Varables {#subsec:cont_ensemble}
--------------------------

We can also apply Theorem 5 to obtain MI estimators that achieve the parametric rate for the case when $\mathbf{X}$ and $\mathbf{Y}$ are continuous. For general $g$, (4) in the main paper indicates that we need $h_{X}^{d_{X}}h_{Y}^{d_{Y}}\propto N^{-1/2}$ for the $O(1/(Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}))$ terms to be $O(1/\sqrt{N})$. We consider the more general case where the parameters may differ for $h_{X}$ and $h_{Y}$. Let $\mathcal{L}_{X}$ and $\mathcal{L}_{Y}$ be sets of real, positive numbers with $|\mathcal{L}_{X}|=L_{X}$ and $|\mathcal{L}_{Y}|=L_{Y}$. For each estimator in the ensemble, choose $l_{X}\in\mathcal{L_{X}}$ and $l_{Y}\in\mathcal{L}_{Y}$ and set $h_{X}(l_{X})=l_{X}N^{-1/(2(d_{X}+d_{Y}))}$ and $h_{Y}(l_{Y})=l_{Y}N^{-1/(2(d_{X}+d_{Y}))}$. Define the matrix $w$ s.t. $\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})=1$. From Theorems 1 and 2, conditions $\mathcal{C}.1$ and $\mathcal{C}.2$ are satisfied if $s\geq d_{X}+d_{Y}$ with $\psi_{i,j}(l_{X},l_{Y})=l_{X}^{i}l_{Y}^{j}$ and $\phi_{i,j}(N)=N^{-(i+j)/(2(d_{X}+d_{Y}))}$ for $0\leq i,j\leq d_{X}+d_{Y}$ s.t. $0<i+j\leq d_{X}+d_{Y}$. The optimal weight $w_{0}$ is calculated using (14) in the main paper. The resulting estimator $$\g{w_{0},1}^{cont}=\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w_{0}(l_{X},l_{Y})\g{h_{X}(l_{X}),h_{Y}(l_{Y})}$$ achieves the parametric MSE rate when $s\geq d_{X}+d_{Y}$.

Again, if the mixed derivatives of the functional $g$ also satisfy the extra condition required for (5) in the main paper, then we can define an estimator that achieves the parametric MSE rate under less strict smoothness assumptions. See Appendix \[subsec:cont\_ensemble2\].

The ODin2 Estimators {#subsec:Odin2}
--------------------

The estimators $\g{w_{0},1}$ and $\g{w_{0},1}^{cont}$ are analogous to the ODin1 estimators in [@moon2016arxiv; @moon2016isit]. In this section, we derive ensemble estimators of MI that achieve the parametric rate under less strict smoothness assumptions on the densities. These estimators are analogous to the ODin2 estimators in [@moon2016arxiv; @moon2016isit].

### Mixed Random Variables {#mixed-random-variables}

We first consider the case where $\mathbf{X}$ is continuous and $\mathbf{Y}$ is discrete. Recall that if $\mathbf{h}_{X|y}=l\mathbf{N}_{y}^{-\beta}$ with $0<\beta<\frac{1}{d_{X}}$ and $l$ a positive number, and if $g\left(t_{1},t_{2}\right)$ has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$ for each $1\leq j,l\leq\lambda$, then the bias of the plug-in estimator for this case is $$\begin{aligned}
\bias\left[\g{h_{X},h_{X|Y}}\right] & =\sum_{\substack{m,n=0\\
i+j+m+n\neq0
}
}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i,j=0}^{r}c_{14,j,i,m,n}\frac{h_{X}^{i}l^{j}N^{-j\beta}}{\left(Nh_{X}^{d_{X}}\right)^{m}\left(l^{d_{X}}N^{1-\beta d_{X}}\right)^{n}}\nonumber \\
 & +O\left(h_{X}^{s}+N^{-s\beta}+\frac{1}{\left(Nh_{X}^{d_{X}}\right)^{\lambda/2}}+\frac{1}{\left(N^{1-\beta d_{X}}\right)^{\lambda/2}}\right).\label{eq:bias_mixed2-1}\end{aligned}$$

Choose $\mathcal{L}$ to be a set of real positive numbers and let $\delta>0$. For each estimator in the ensemble, set $h_{X}(l)=lN^{-1/(d_{X}+\delta)}$ and $\mathbf{h}_{X|y}(l)=l\mathbf{N}_{y}^{-1/(d_{X}+\delta)}$ where $l\in\mathcal{L}$. This ensures that the final terms in (\[eq:bias\_mixed2-1\]) are $O(1/\sqrt{N})$ if $s\geq(d_{X}+\delta)/2$ and $\lambda\geq d_{X}/\delta+1$. Define $\g{w,2}$ as in (12) in the main paper with the chosen values of $h_{X}(l)$ and $\mathbf{h}_{X|y}(l)$. Theorem 5 can be applied in this case as conditions $\mathcal{C}.1$ and $\mathcal{C}.2$ are satisfied with $\phi_{i,m}(N)=N^{\frac{-i-m\delta}{d_{X}+\delta}}$ and $\psi_{i,m}(l)=l^{i-md_{X}}$ for $i\in\{0,\dots,r\}$, $m\in\{0,\dots\left\lfloor \lambda/2\right\rfloor \}$, and $0\leq i+m\delta\leq(d_{X}+\delta)/2$. Then if $s\geq(d_{X}+\delta)/2$ and $\lambda\geq d_{X}/\delta+1$, the MSE of the optimally weighted estimator $\g{w_{0},2}$ is $O(1/N)$. Then since $\delta$ can be chosen arbitrarily close to zero, the parametric rate can be achieved theoretically as long as $s>d_{X}/2$.

The analogous divergence functional estimators for $\g{w_{0},1}$ and $\g{w_{0},2}$ in [@moon2016arxiv; @moon2016isit] were referred to as the ODin1 and ODin2 estimators, respectively, where ODin stands for **O**ptimally Weighted **Di**stributional Fu**n**ctional estimators. The ODin2 estimator has better statistical properties as the parametric rate is guaranteed under less restrictive smoothness assumptions on the densities. On the other hand, the number of parameters required for the optimization problem in (14) in the main paper is larger for the ODin2 estimator than the ODin1 estimator. In theory, this could lead to larger variance although this wasn’t necessarily true in practice according to the experiments in [@moon2016arxiv].

### Continuous Random Variables {#subsec:cont_ensemble2}

We now consider the case where both $\mathbf{X}$ and $\mathbf{Y}$ are continuous. Again, if $g\left(t_{1},t_{2}\right)$ has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$ for each $1\leq k,l,\leq\lambda$, then the bias of $\gt$ is $$\begin{aligned}
\bias\left[\gt\right] & = & \sum_{\substack{m,n=0\\
i+j+m+n\neq0
}
}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i,j=0}^{r}c_{11,j,i,m,n}\frac{h_{X}^{i}h_{Y}^{j}}{\left(Nh_{X}^{d_{X}}\right)^{m}\left(Nh_{Y}^{d_{Y}}\right)^{n}}\nonumber \\
 &  & +\sum_{m=1}^{\left\lfloor \lambda/2\right\rfloor }\sum_{i=0}^{r}\sum_{j=0}^{r}c_{13,m,n,j}h_{X}^{i}h_{Y}^{j}/\left(Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\right)^{m}\nonumber \\
 &  & +O\left(h_{X}^{s}+h_{Y}^{s}+1/\left(Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\right)^{\lambda/2}\right).\label{eq:bias2-1}\end{aligned}$$

Set $\delta>0$ and choose $h_{X}(l_{X})=l_{X}N^{-1/(d_{X}+d_{Y}+\delta)}$ and $h_{Y}(l_{Y})=l_{Y}N^{-1/(d_{X}+d_{Y}+\delta)}$. Then conditions $\mathcal{C}.1$ and $\mathcal{C}.2$ are satisfied if $s\geq(d_{X}+d_{Y}+\delta)/2$ and $\lambda\geq(d_{X}+d_{Y}+\delta)/\delta$ with $\psi_{1,i,j,m,n}(l_{X},l_{Y})=l_{X}^{i-md_{X}}l_{Y}^{j-nd_{Y}}$ and $\phi_{1,i,j,m,n}(N)=N^{-\frac{i+j+m(d_{Y}+\delta)+n(d_{X}+\delta)}{d_{X}+d_{Y}+\delta}}$ for $0<i+j+m(d_{Y}+\delta)+n(d_{X}+\delta)\leq\frac{d_{X}+d_{Y}+\delta}{2}$ and the terms $\psi_{2,i,j,m}(l_{X},l_{Y})=l_{X}^{i-md_{X}}l_{Y}^{j-md_{Y}}$ and $\phi_{2,i,j,m}(N)=N^{-\frac{i+j+m\delta}{d_{X}+d_{Y}+\delta}}$ for $m\geq1$ and $i+j+m\delta\leq\frac{d_{X}+d_{Y}+\delta}{2}.$ The optimal weight $w_{0}$ is again calculated using (14) in the main paper and the resulting estimator $\g{w_{0},2}^{cont}$ achieves the parametric MSE convergence rate when $s\geq(d_{X}+d_{Y}+\delta)/2$. Since $\delta$ can be chosen arbitrarily close to zero, the parametric rate can be achieved theoretically as long as $s>(d_{X}+d_{Y})/2$. $\g{w_{0},2}^{cont}$ is the ODin2 estimator for continuous random variables.

Proof of Theorem 1 (Bias) {#sec:biasProof}
=========================

The proof of the bias results in Theorem 1 share some similarities with the proof of the bias results for the divergence functional estimators in in [@moon2016arxiv]. The primary differences deal with the product of the marginal KDEs that appear in the expansion of the bias terms.

The bias of $\gt$ can be expressed as $$\begin{aligned}
\bias\left[\gt\right] & = & \bE\left[g\left(\frac{\ft X(\mathbf{X})\ft Y(\mathbf{Y})}{\ft Z(\mathbf{X},\mathbf{Y})}\right)-g\left(\frac{f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})}{f_{XY}(\mathbf{X},\mathbf{Y})}\right)\right]\nonumber \\
 & = & \bE\left[g\left(\frac{\ft X(\mathbf{X})\ft Y(\mathbf{Y})}{\ft Z(\mathbf{X},\mathbf{Y})}\right)-g\left(\frac{\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]}{\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})}\right)\right]\nonumber \\
 &  & +\bE\left[g\left(\frac{\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]}{\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})}\right)-g\left(\frac{f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})}{f_{XY}(\mathbf{X},\mathbf{Y})}\right)\right],\label{eq:gsplit}\end{aligned}$$ where $\mathbf{X}$ and $\mathbf{Y}$ are drawn jointly from $f_{XY}$. We can view these terms as a variance-like component (the first term) and a bias-like component, where the respective Taylor series expansions depend on variance-like or bias-like terms of the KDEs.

We first consider the bias-like term, i.e. the second term in (\[eq:gsplit\]). The Taylor series expansion of $g\left(\frac{\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]}{\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})}\right)$ around $f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})$ and $f_{XY}(\mathbf{X},\mathbf{Y})$ gives an expansion with terms of the form of $$\begin{aligned}
\bias_{\mathbf{Z}}^{i}\left[\ft X(\mathbf{X})\ft Y(\mathbf{Y})\right] & = & \left(\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]-f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})\right)^{i},\nonumber \\
\bias_{\mathbf{Z}}^{i}\left[\ft Z(\mathbf{X},\mathbf{Y})\right] & = & \left(\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})-f_{XY}(\mathbf{X},\mathbf{Y})\right)^{i}.\label{eq:biasterms}\end{aligned}$$

Since we are not doing boundary correction, we need to consider separately the cases when $\mathbf{Z}$ is in the interior of the support $\mathcal{S}_{X}\times\mathcal{S}_{Y}$ and when $\mathbf{Z}$ is close to the boundary of the support. For precise definitions, a point $Z=(X,Y)\in\mathcal{S}_{X}\times\mathcal{S}_{Y}$ is in the interior of $\mathcal{S}_{X}\times\mathcal{S}_{Y}$ if for all $Z^{'}\notin\mathcal{S}_{X}\times\mathcal{S}_{Y}$, $K_{X}\left(\frac{X-X^{'}}{h_{X}}\right)K_{Y}\left(\frac{Y-Y^{'}}{h_{Y}}\right)=0$, and a point $Z\in\mathcal{S}_{X}\times\mathcal{S}_{Y}$ is near the boundary of the support if it is not in the interior.

It can be shown by Taylor series expansions of the probability densities that for $\mathbf{Z}=(\mathbf{X},\mathbf{Y})$ drawn from $f_{XY}$ in the interior of $\mathcal{S}_{X}\times\mathcal{S}_{Y}$, then $$\begin{aligned}
\ez X\left[\ft X(\mathbf{X})\right] & = & f_{X}(\mathbf{X})+\sum_{j=1}^{\left\lfloor s/2\right\rfloor }c_{X,j}(\mathbf{X})h_{X}^{2j}+O\left(h_{X}^{s}\right),\label{eq:E_fx}\\
\ez Y\left[\ft Y(\mathbf{Y})\right] & = & f_{Y}(\mathbf{Y})+\sum_{j=1}^{\left\lfloor s/2\right\rfloor }c_{Y,j}(\mathbf{Y})h_{Y}^{2j}+O\left(h_{Y}^{s}\right),\nonumber \\
\ez{X,Y}\left[\ft Z(\mathbf{Z})\right] & = & f_{XY}(\mathbf{X},\mathbf{Y})+\sum_{\substack{i=0\\
i+j\neq0
}
}^{\left\lfloor s/2\right\rfloor }\sum_{j=0}^{\left\lfloor s/2\right\rfloor }c_{XY,i,j}(\mathbf{X},\mathbf{Y})h_{X}^{2i}h_{Y}^{2j}+O\left(h_{X}^{s}+h_{Y}^{s}\right).\nonumber \end{aligned}$$

For a point near the boundary of the support, we extend the expectation beyond the support of the density. As an example if $\mathbf{X}$ is near the boundary of $\mathcal{S}_{X}$, then we get $$\begin{aligned}
\mathbb{E}_{\mathbf{X}}\left[\ft i(\mathbf{X})\right]-f_{i}(\mathbf{X}) & = & \frac{1}{h_{X}^{d_{X}}}\int_{V:V\in\mathcal{S}_{X}}K_{X}\left(\frac{\mathbf{X}-V}{h_{X}}\right)f_{X}(V)dV-f_{X}(\mathbf{X})\nonumber \\
 & = & \left[\frac{1}{h_{X}^{d_{X}}}\int_{V:K_{X}\left(\frac{\mathbf{X}-V}{h_{X}}\right)>0}K_{X}\left(\frac{\mathbf{X}-V}{h_{X}}\right)f_{X}(V)dV-f_{X}(\mathbf{X})\right]\nonumber \\
 &  & -\left[\frac{1}{h_{X}^{d_{X}}}\int_{V:V\notin\mathcal{S}_{X}}K_{X}\left(\frac{\mathbf{X}-V}{h_{X}}\right)f_{X}(V)dV\right]\nonumber \\
 & = & T_{1,X}(\mathbf{X})-T_{2,X}(\mathbf{X}).\label{eq:Tdiff}\end{aligned}$$ We only evaulate the density $f_{X}$ and its derivatives at points within the support when we take its Taylor series expansion. Thus the exact manner in which we define the extension of $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}(\mathbf{X})$ gives an expression of the form of (\[eq:E\_fx\]). For the $T_{2,X}(\mathbf{X})$ term, we can use multi-index notation on the expansion of $f_{X}$ to show that $$\begin{aligned}
T_{2,X}(\mathbf{X}) & =\left[\frac{1}{h_{X}^{d_{X}}}\int_{V:V\notin\mathcal{S}_{X}}K_{X}\left(\frac{\mathbf{X}-V}{h_{X}}\right)f_{X}(V)dV\right]\\
 & =\int_{u:h_{X}u+\mathbf{X}\notin\mathcal{S}_{X},K_{X}(u)>0}K_{X}(u)f_{X}(\mathbf{X}+h_{X}u)du\\
 & =\sum_{|\alpha|\leq r}\frac{h_{X}^{|\alpha|}}{\alpha!}\int_{u:h_{X}u+\mathbf{X}\notin\mathcal{S}_{X},K_{X}(u)>0}K_{X}(u)D^{\alpha}f_{X}(\mathbf{X})u^{\alpha}du+o(h_{X}^{r}).\end{aligned}$$ Then since the $|\alpha|$th derivative of $f_{X}$ is $r-|\alpha|$ times differentiable, we apply the condition in assumption $\mathcal{A}.5$ to obtain $$\bE\left[T_{2,X}(\mathbf{X})\right]=\sum_{i=1}^{r}e_{i}h_{X}^{i}+o\left(h_{X}^{r}\right).$$ Similar expressions can be found for $\ft Y$ and $\ft Z$ and for when (\[eq:Tdiff\]) is raised to a power $t$. Applying this result gives for the second term in (\[eq:gsplit\]), $$\sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{10,i,j}h_{X}^{i}h_{Y}^{j}+O\left(h_{X}^{s}+h_{Y}^{s}\right).\label{eq:varterm_bound}$$

For the first term in (\[eq:gsplit\]), a Taylor series expansion of $g\left(\frac{\ft X(\mathbf{X})\ft Y(\mathbf{Y})}{\ft Z(\mathbf{X},\mathbf{Y})}\right)$ around $\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]$ and $\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})$ gives an expansion with terms of the form of $$\begin{aligned}
\et Z^{q}(\mathbf{Z}) & = & \left(\ft Z(\mathbf{Z})-\ez Z\left[\ft Z(\mathbf{Z})\right]\right)^{q},\nonumber \\
\ett XY^{q}(\mathbf{Z}) & = & \left(\ft X(\mathbf{X})\ft Y(\mathbf{Y})-\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]\right)^{q}.\label{eq:varterms}\end{aligned}$$ We can take the expected value of these expressions to obtain terms of the form of $$\frac{1}{Nh_{X}^{d_{X}}},\,\frac{1}{Nh_{Y}^{d_{Y}}},\,\frac{1}{N^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}},\,\frac{1}{Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}}\label{eq:variance_terms}$$ and their respective powers. This can be seen for $\ett XY^{q}(\mathbf{Z})$ as follows. Define $$\begin{aligned}
\mathbf{V}_{i,j}(\mathbf{Z}) & =K_{X}\left(\frac{\mathbf{X}_{i}-\X}{h_{X}}\right)K_{Y}\left(\frac{\Y_{j}-\Y}{h_{Y}}\right)-\ez X\left[K_{X}\left(\frac{\mathbf{X}_{i}-\X}{h_{X}}\right)\right]\ez Y\left[K_{Y}\left(\frac{\Y_{j}-\Y}{h_{Y}}\right)\right]\\
 & =\eta_{ij}(\Z)-\ez X\left[\eta_{i}(\X)\right]\ez Y\left[\eta_{j}^{'}(\Y)\right].\end{aligned}$$ We can then write $$\ett XY(\mathbf{Z})=\frac{1}{N^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathbf{V}_{i,j}(\Z).$$ The binomial theorem then gives $$\ez Z\left[\mathbf{V}_{i,j}^{k}(\Z)\right]=\sum_{l=0}^{k}\binom{k}{l}\ez Z\left[\eta_{ij}^{l}(\Z)\right]\left(\ez X\left[\eta_{i}(\X)\right]\ez Y\left[\eta_{j}^{'}(\Y)\right]\right)^{k-l}.\label{eq:exp_V}$$ By using a similar Taylor series analysis as before, for $\Z$ in the interior, $$\ez Z\left[\eta_{ij}^{l}(\Z)\right]=h_{X}^{d_{X}}h_{Y}^{d_{Y}}\sum_{m,n=0}^{\left\lfloor s/2\right\rfloor }c_{XY,2,m,n,l}(\Z)h_{X}^{2m}h_{Y}^{2n}+O\left(h_{X}^{2d_{X}}h_{Y}^{d_{Y}}+h_{X}^{d_{X}}h_{Y}^{2d_{Y}}\right).$$ Combining this with (\[eq:E\_fx\]) and (\[eq:exp\_V\]) gives $$\ez Z\left[\mathbf{V}_{i,j}^{k}(\Z)\right]=h_{X}^{d_{X}}h_{Y}^{d_{Y}}\sum_{m,n=0}^{\left\lfloor s/2\right\rfloor }c_{XY,3,m,n,k}(\X)h_{X}^{2m}h_{Y}^{2n}+O\left(h_{X}^{2d_{X}}h_{Y}^{d_{Y}}+h_{X}^{d_{X}}h_{Y}^{2d_{Y}}\right),\label{eq:exp_Vk}$$ where the constants depend on the densities, their derivatives, and the moments of the kernels. As an example, let $q=2$. Then due to the independence between the $\Z_{i}$ samples, $$\begin{aligned}
\ez Z\left[\ett XY^{2}(\mathbf{Z})\right] & =\frac{1}{N^{4}h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}\sum_{i,j,m,n=1}^{N}\ez Z\left[\mathbf{V}_{i,j}(\Z)\mathbf{V}_{m,n}(\Z)\right]\\
 & =\frac{1}{N^{2}h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}\ez Z\left[\mathbf{V}_{i,j}^{2}(\Z)\right]+\frac{(N-1)}{N^{2}h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}\ez Z\left[\mathbf{V}_{i,j}(\Z)\mathbf{V}_{i,n}(\Z)\right]\\
 & =\frac{1}{N^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}}\sum_{m,n=0}^{\left\lfloor s/2\right\rfloor }c_{XY,3,m,n,2}(\X)h_{X}^{2m}h_{Y}^{2n}+\sum_{m,n=0}^{\left\lfloor s/2\right\rfloor }\sum_{\substack{i,j=0\\
i+j\neq0
}
}^{1}c_{XY,4,m,n,i,j}(\X)\frac{h_{X}^{2m}h_{Y}^{2n}}{Nh_{X}^{id_{X}}h_{Y}^{jd_{Y}}}+O\left(\frac{1}{N}\right),\end{aligned}$$ where the last step follows from (\[eq:exp\_Vk\]) and a similar analysis of $\ez Z\left[\mathbf{V}_{i,j}(\Z)\mathbf{V}_{i,n}(\Z)\right]$. For $q>2$, it can be shown that if $n(q)$ is the set of integer divisors of $q$ including 1 but excluding $q$, then $$\ez Z\left[\ett XY^{q}(\mathbf{Z})\right]=\sum_{i,j=0}^{\left\lfloor s/2\right\rfloor }\left(\sum_{n\in n(q)}\frac{c_{XY,5,i,j,q,n}(\Z)}{\left(N^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}\right)^{q-n}}+\sum_{\substack{m\in n(q)\cup\{q\}\\
n\in n(q)\cup\{q\}\\
m+n\neq2q
}
}\frac{c_{XY,6,i,j,q,m,n}(\Z)}{\left(Nh_{X}^{d_{X}}\right)^{q-n}\left(Nh_{Y}^{d_{Y}}\right)^{q-m}}\right)h_{X}^{2i}h_{Y}^{2j}+O\left(\frac{1}{N}\right).$$ A similar procedure can be used to find the expression for $\ez Z\left[\et Z^{q}(\mathbf{Z})\right]$. When $\Z$ is near the boundary of the supposrt, we can obtain similar expressions by following a similar procedure as in the derivation of (\[eq:varterm\_bound\]). This results in powers of $h_{X}^{m}h_{Y}^{n}$ instead of $h_{X}^{2m}h_{Y}^{2n}$.

For general functionals $g$, we can only guarantee that the mixed derivatives of $g$ evaluated at $\ez X\left[\ft X(\mathbf{X})\right]\ez Y\left[\ft Y(\mathbf{Y})\right]$ and $\ez{X,Y}\ft Z(\mathbf{X},\mathbf{Y})$ converge to the mixed derivative evaluated at $f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})$ and $f_{XY}(\mathbf{X},\mathbf{Y})$ at some rate $o(1)$. Thus we are left with the following terms in the bias: $$o\left(\frac{1}{Nh_{X}^{d_{X}}}+\frac{1}{Nh_{Y}^{d_{Y}}}\right)$$ However, if we know that $g\left(t_{1},t_{2}\right)$ has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$, then by the generalized binomial theorem, we find that $$\left(\ez{\mathbf{X}}\ft X(\mathbf{X})\right)^{\alpha}=\sum_{m=0}^{\infty}\binom{\alpha}{m}f_{X}^{\alpha-m}(\mathbf{X})\left(\sum_{j=1}^{\left\lfloor s/2\right\rfloor }c_{i,j}(\mathbf{X})h_{X}^{2j}+O\left(h_{X}^{s}\right)\right)^{m}.$$ A similar result holds for $\left(\ez{\mathbf{Y}}\ft Y(\mathbf{Y})\right)^{\alpha}$ and $\left(\ez Z\ft Z(\mathbf{Z})\right)^{\alpha}$. Combining these expressions with \[eq:variance\_terms\] completes the proof.

Proof of Theorem 2 (Variance) {#sec:VarProof}
=============================

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

(Efron-Stein Inequality) Let $\mathbf{X}_{1},\dots,\mathbf{X}_{n},\mathbf{X}_{1}^{'},\dots,\mathbf{X}_{n}^{'}$ be independent random variables on the space $\mathcal{S}$. Then if $f:\mathcal{S}\times\dots\times\mathcal{S}\rightarrow\mathbb{R}$, we have that $$\var\left[f(\mathbf{X}_{1},\dots,\mathbf{X}_{n})\right]\leq\frac{1}{2}\sum_{i=1}^{n}\bE\left[\left(f(\mathbf{X}_{1},\dots,\mathbf{X}_{n})-f(\mathbf{X}_{1},\dots,\mathbf{X}_{i}^{'},\dots,\mathbf{X}_{n})\right)^{2}\right].$$

In this case we consider the samples $\left\{ \mathbf{Z}_{1},\dots,\mathbf{Z}_{N}\right\} $ and $\left\{ \mathbf{Z}_{1}^{'},\Z_{2}\dots,\mathbf{Z}_{N}\right\} $ and the respective estimators $\gt$ and $\gt^{'}$. By the triangle inequality, $$\begin{aligned}
\left|\gt-\gt^{'}\right| & \leq & \frac{1}{N}\left|g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)\right|\nonumber \\
 &  & +\frac{1}{N}\sum_{j=2}^{N_{2}}\left|g\left(\frac{\ft X(\mathbf{X}_{j})\ft Y(\mathbf{Y}_{j})}{\ft Z(\mathbf{X}_{j},\mathbf{Y}_{j})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{j})\ft Y^{'}(\mathbf{Y}_{1})}{\ft Z^{'}(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right|.\label{eq:triangle}\end{aligned}$$ By the Lipschitz condition on $g$, the first term in (\[eq:triangle\]) can be decomposed into terms of the form of $$\left|\ft Z(\mathbf{Z}_{1})-\ft Z(\mathbf{Z}_{1}^{'})\right|,$$ $$\left|\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})-\ft X(\mathbf{X}_{1}^{'})\ft Y^{'}(\mathbf{Y}_{1})\right|.$$ By making a substitution in the expectation, it can be shown that $$\bE\left[\left|\ft Z(\mathbf{Z}_{1})-\ft Z(\mathbf{Z}_{1}^{'})\right|^{2}\right]\leq2||K_{X}\cdot K_{Y}||_{\infty}^{2}.$$ For the product of the marginal KDEs, we have that $$\begin{aligned}
\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1}) & = & \frac{1}{M^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}}\sum_{i=2}^{N}\sum_{j=2}^{N}K_{X}\left(\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}-\mathbf{Y}_{j}}{h_{Y}}\right)\\
 & = & \frac{1}{M}\ft Z(\mathbf{Z}_{1})+\frac{1}{M^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}}\sum_{i\neq j}K_{X}\left(\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}-\mathbf{Y}_{j}}{h_{Y}}\right).\end{aligned}$$ By applying the triangle inequality, Jensen’s inequality, and similar substitutions, we get $$\begin{aligned}
\bE\left[\left|\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})-\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})\right|^{2}\right] & \leq & \bE\left[\frac{2}{M^{2}}\left|\ft Z(\mathbf{Z}_{1})-\ft Z(\mathbf{Z}_{1}^{'})\right|^{2}\right]\\
 &  & +\frac{2(M-1)}{M^{3}h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}\times\\
 &  & \sum_{i\neq j}\bE\left[\left(K_{X}\left(\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}-\mathbf{Y}_{j}}{h_{Y}}\right)\right.\right.\\
 &  & \left.\left.-K_{X}\left(\frac{\mathbf{X}_{1}^{'}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}^{'}-\mathbf{Y}_{j}}{h_{Y}}\right)\right)^{2}\right]\\
 & \leq & \frac{4+2(M-1)^{2}}{M^{2}}||K_{X}\cdot K_{Y}||^{2}.\end{aligned}$$

For the second term in (\[eq:triangle\]), it can be shown that $$\begin{aligned}
\bE\left[\left|\ft Z(\mathbf{Z}_{i})-\ft Z^{'}(\mathbf{Z}_{i})\right|^{2}\right] & = & \frac{1}{M^{2}h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}\bE\left[\left(K_{X}\left(\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}-\mathbf{Y}_{j}}{h_{Y}}\right)\right.\right.\\
 &  & \left.\left.-K_{X}\left(\frac{\mathbf{X}_{1}^{'}-\mathbf{X}_{i}}{h_{X}}\right)K_{Y}\left(\frac{\mathbf{Y}_{1}^{'}-\mathbf{Y}_{j}}{h_{Y}}\right)\right)^{2}\right]\\
 & \leq & \frac{2||K_{X}\cdot K_{Y}||_{\infty}^{2}}{M^{2}}.\end{aligned}$$ By a similar approach, $$\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})-\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})$$ $$\begin{aligned}
 & = & \ft Z(\mathbf{Z}_{i})-\ft Z^{'}(\mathbf{Z}_{i})+\frac{1}{M^{2}h_{X}^{d_{X}}h_{Y}^{d_{Y}}}\left(\sum_{\substack{n=2\\
n\neq i
}
}K_{Y}\left(\frac{\mathbf{Y}_{i}-\mathbf{Y}_{n}}{h_{Y}}\right)\left(K_{X}\left(\frac{\mathbf{X}_{i}-\mathbf{X}_{1}}{h_{X}}\right)-K_{X}\left(\frac{\mathbf{X}_{i}-\mathbf{X}_{1}^{'}}{h_{X}}\right)\right)\right.\\
 &  & \left.+\sum_{\substack{n=2\\
n\neq i
}
}K_{X}\left(\frac{\mathbf{X}_{i}-\mathbf{X}_{n}}{h_{X}}\right)\left(K_{Y}\left(\frac{\mathbf{Y}_{i}-\mathbf{Y}_{1}}{h_{Y}}\right)-K_{Y}\left(\frac{\mathbf{Y}_{i}-\mathbf{Y}_{1}^{'}}{h_{Y}}\right)\right)\right),\end{aligned}$$ $$\implies\bE\left[\left|\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})-\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})\right|^{2}\right]\leq6||K_{X}\cdot K_{Y}||_{\infty}^{2}\left(\frac{1}{M^{2}}+\frac{(M-2)^{2}}{M^{4}}\right)$$

We can then apply the Cauchy Schwarz inequality to bound the square of the second term in (\[eq:triangle\]) to get $$\bE\left[\left(\sum_{j=2}^{N_{2}}\left|g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{1})\ft Y^{'}(\mathbf{Y}_{1})}{\ft Z^{'}(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right|\right)^{2}\right]\leq14C_{g}^{2}||K_{X}\cdot K_{Y}||_{\infty}^{2}.$$ Applying Jensen’s inequality in conjunction with these results gives $$\bE\left[\left|\gt-\gt^{'}\right|^{2}\right]\leq\frac{44C_{g}^{2}||K_{X}\cdot K_{Y}||_{\infty}^{2}}{N^{2}}.$$ Applying the Efron-Stein inequality finishes the proof.

Theory for Mixed Random Variables\[sec:MixedProofs\]
====================================================

Proof of Theorem 3 (Bias)
-------------------------

Let $\mathbf{h}_{X|y}=l\mathbf{N}_{y}^{-\beta}$ for some positive $l$ and $0<\beta<\frac{1}{d_{X}}$. Under assumptions $\mathcal{A}.0-\mathcal{A}.5$, we prove that for general $g$, the bias of the plug-in estimator $\g{h_{X},h_{X|Y}}$ $$\begin{aligned}
\bias\left[\g{h_{X},h_{X|Y}}\right] & = & \sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{13,i,j}h_{X}^{i}l^{j}N^{-j\beta}+\frac{c_{14,X}}{Nh_{X}^{d_{X}}}+\frac{c_{14,y}}{l^{d_{X}}N^{1-\beta d_{X}}}\nonumber \\
 &  & +O\left(h_{X}^{s}+N^{-s\beta}+\frac{1}{Nh_{X}^{d_{X}}}+\frac{1}{N^{1-\beta d_{X}}}+\frac{1}{N}\right).\label{eq:bias_mixed1-1}\end{aligned}$$ Furthermore, if $g\left(t_{1},t_{2}\right)$ has $j,l$-th order mixed derivatives $\frac{\partial^{j+l}}{\partial t_{1}^{j}\partial t_{2}^{l}}$ that depend on $t_{1}$ and $t_{2}$ only through $t_{1}^{\alpha}t_{2}^{\beta}$ for some $\alpha,\beta\in\mathbb{R}$, then for any positive integer $\lambda\geq2$, the bias is $$\begin{aligned}
\bias\left[\g{h_{X},h_{X|Y}}\right] & = & \sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{13,i,j}h_{X}^{i}l^{j}N^{-j\beta}+\sum_{j=1}^{\lambda/2}\sum_{i=1}^{\lambda/2}\sum_{m=0}^{r}\sum_{n=0}^{r}c_{14,j,i,m,n}\frac{h_{X}^{m}l^{n}N^{-n\beta}}{\left(Nh_{X}^{d_{X}}\right)^{j}\left(l^{d_{X}}N^{1-\beta d_{X}}\right)^{i}}\nonumber \\
 &  & +\sum_{j=1}^{\lambda/2}\sum_{m=0}^{r}\sum_{n=0}^{r}\left(c_{14,m,n,j,X}\frac{h_{X}^{m}l^{n}N^{-n\beta}}{\left(Nh_{X}^{d_{X}}\right)^{j}}+c_{14,m,n,j,Y}\frac{h_{X}^{m}l^{n}N^{-n\beta}}{\left(l^{d_{X}}N^{1-\beta d_{X}}\right)^{j}}\right)\nonumber \\
 &  & +O\left(h_{X}^{s}+N^{-s\beta}+\frac{1}{\left(Nh_{X}^{d_{X}}\right)^{\lambda/2}}+\frac{1}{\left(N^{1-\beta d_{X}}\right)^{\lambda/2}}+\frac{1}{N}\right).\label{eq:bias_mixed2-1}\end{aligned}$$

We only prove (\[eq:bias\_mixed1-1\]) as the proof of (\[eq:bias\_mixed2-1\]) is identical. The bias of $\g{h_{X},h_{X|Y}}$ is $$\begin{aligned}
\bias\left[\g{h_{X},h_{X|Y}}\right] & = & \bE\left[\g{h_{X},h_{X|Y}}\right]-G(\mathbf{X};\mathbf{Y})\\
 & = & \bE\left[\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\g{h_{X},h_{X|y}}-g\left(\frac{f_{X}(\mathbf{X})}{f_{X|Y}(\mathbf{X}|\mathbf{Y})}\right)\right]\\
 & = & \bE\left[\bE\left[\left.\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\g{h_{X},h_{X|y}}-g\left(\frac{f_{X}(\mathbf{X})}{f_{X|Y}(\mathbf{X}|\mathbf{Y})}\right)\right|\mathbf{Y},\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right]\\
 & = & \bE\left[\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\bE\left[\left.\left(\g{h_{X},h_{X|y}}-g\left(\frac{f_{X}(\mathbf{X})}{f_{X|Y}(\mathbf{X}|\mathbf{Y})}\right)\right)\right|\mathbf{Y},\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right]\\
 & = & \bE\left[\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}\bias\left[\left.\g{h_{X},h_{X|y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right],\end{aligned}$$ where we use the law of total expectation and the fact that $\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}=1$. Let $\mathbf{h}_{X|y}=l\mathbf{N}_{y}^{-\beta}$ for some positive $l$ and $0<\beta<\frac{1}{d_{X}}$. From Theorem 1, the conditional bias of $\g{h_{X},h_{X|y}}$ given $\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}$ is $$\begin{aligned}
\bias\left[\left.\g{h_{X},h_{X|y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right] & = & \sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{10,i,j}h_{X}^{i}\mathbf{h}_{X|y}^{j}+\frac{c_{11,X}}{\mathbf{N}_{y}h_{X}^{d_{X}}}+\frac{c_{11,y}}{\mathbf{N}_{y}\mathbf{h}_{X|y}^{d_{X}}}\nonumber \\
 &  & +O\left(h_{X}^{s}+\mathbf{h}_{X|y}^{s}+\frac{1}{\mathbf{N}_{y}h_{X}^{d_{X}}}+\frac{1}{\mathbf{N}_{y}\mathbf{h}_{X|y}^{d_{X}}}\right)\nonumber \\
 &  & =\sum_{\substack{j=0\\
i+j\neq0
}
}^{r}\sum_{i=0}^{r}c_{10,i,j}h_{X}^{i}l^{j}\mathbf{N}_{y}^{-j\beta}+\frac{c_{11,X}}{\mathbf{N}_{y}h_{X}^{d_{X}}}+\frac{c_{11,y}}{l^{d_{X}}\mathbf{N}_{y}^{1-\beta d_{X}}}\nonumber \\
 &  & +O\left(h_{X}^{s}+\mathbf{N}_{y}^{-s\beta}+\frac{1}{\mathbf{N}_{y}h_{X}^{d_{X}}}+\frac{1}{\mathbf{N}_{y}^{1-\beta d_{X}}}\right).\label{eq:bias_cond}\end{aligned}$$ $\mathbf{N}_{y}$ is a binomial random variable Multiplying (\[eq:bias\_cond\]) by $\mathbf{N}_{y}$ results in terms of the form of $\mathbf{N}_{y}^{1-\gamma}$ with $\gamma\geq0$. $\mathbf{N}_{y}$ is a binomial random variable with parameter $f_{Y}(y)$,$N$ trials, and mean $Nf_{Y}(y)$. We can compute the fractional moments of a binomial random variable by using the generalized binomial theorem to obtain (see the main paper) $$\begin{aligned}
\bE\left[\mathbf{N}_{y}^{\alpha}\right] & = & \sum_{i=0}^{\infty}\left(\begin{array}{c}
\alpha\\
i
\end{array}\right)\left(Nf_{Y}(y)\right)^{\alpha-i}\bE\left[\left(\mathbf{N}_{Y}-Nf_{Y}(y)\right)^{i}\right]\\
 & = & \sum_{i=0}^{\infty}\left(\begin{array}{c}
\alpha\\
i
\end{array}\right)\left(Nf_{Y}(y)\right)^{\alpha-i}\sum_{n=0}^{\left\lfloor i/2\right\rfloor }c_{n,i}(f_{Y}(y))N^{n}\\
 & = & \sum_{i=0}^{\infty}\left(\begin{array}{c}
\alpha\\
i
\end{array}\right)f_{Y}(y)^{\alpha-i}\sum_{n=0}^{\left\lfloor i/2\right\rfloor }c_{n,i}(f_{Y}(y))N^{\alpha-i+n},\end{aligned}$$ where we use the following expression for the $i$-th central moment of a binomial random variable derived by Riordan [@riordan1937moment]: $$\bE\left[\left(\mathbf{N}_{Y}-Nf_{Y}(y)\right)^{i}\right]=\sum_{n=0}^{\left\lfloor i/2\right\rfloor }c_{n,i}(f_{Y}(y))N^{n}.$$ If $\alpha=1-\gamma$, then dividing by $N$ results in terms of the form of $N^{-\gamma-i+n}$. Since $n\leq\left\lfloor i/2\right\rfloor $, $-\gamma-i+n$ is always less than zero and is only greater than $-1$ if $i=0$. This completes the proof.

Proof of Theorem 4 (Variance) {#subsec:varproof}
-----------------------------

As for the bias, we assume that $\mathbf{h}_{X|y}=l\mathbf{N}_{y}^{-\beta}$ for some positive $l$ and $0<\beta<\frac{1}{d_{X}}$. By the law of total variance, we have $$\var\left[\g{h_{X},h_{X|Y}}\right]=\bE\left[\var\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right]+\var\left[\bE\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right].\label{eq:total_var}$$ Note that given all of the $\mathbf{Y}_{i}$’s, the estimators $\g{h_{X},h_{X|y}}$ are all independent since they use different sets of $\mathbf{X}_{i}$’s for each $y$. By Theorem 2, we have $$\begin{aligned}
\var\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right] & = & O\left(\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}^{2}}{N^{2}}\cdot\frac{1}{\mathbf{N}_{y}}\right)\\
 & = & O\left(\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N^{2}}\right).\end{aligned}$$ Taking the expectation wrt $\mathbf{Y}_{1},\dots\mathbf{Y}_{N}$ then gives $O\left(\frac{1}{N}\right)$ for the first term in (\[eq:total\_var\]).

For the second term in (\[eq:total\_var\]), from (\[eq:bias\_cond\]) we have that for general $g$ $$\begin{aligned}
\bE\left[\left.\g{h_{X},h_{X|y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right] & = & O\left(\sum_{j=0}^{r}\mathbf{N}_{y}^{-j\beta}+\frac{1}{\mathbf{N}_{y}}+\mathbf{N}_{y}^{-s\beta}+\mathbf{N}_{y}^{1-\beta d_{X}}\right)\\
 & = & O\left(f\left(\mathbf{N}_{y}\right)\right).\end{aligned}$$ By the Efron-Stein inequality, we have that if $\mathbf{N}_{y}^{'}$ is an independent and identically distributed realization of $\mathbf{N}_{y}$, then $$\begin{aligned}
\var\left[\sum_{y\in\mathcal{S}_{Y}}\frac{\mathbf{N}_{y}}{N}f\left(\mathbf{N}_{y}\right)\right] & \leq & \frac{1}{2N^{2}}\sum_{y\in\mathcal{S}_{Y}}\bE\left[\left(\mathbf{N}_{y}f\left(\mathbf{N}_{y}\right)-\mathbf{N}_{y}^{'}f\left(\mathbf{N}_{y}^{'}\right)\right)^{2}\right]\nonumber \\
 & = & O\left(\frac{1}{N^{2}}\bE\left[\left(\mathbf{N}_{y}f\left(\mathbf{N}_{y}\right)-\mathbf{N}_{y}^{'}f\left(\mathbf{N}_{y}^{'}\right)\right)^{2}\right]\right)\nonumber \\
 & = & O\left(\frac{1}{N^{2}}\var\left[\mathbf{N}_{y}f\left(\mathbf{N}_{y}\right)\right]\right),\label{eq:var_cond}\end{aligned}$$ where the second step follows from the fact that $\mathcal{S}_{Y}$ is finite and the last step follows from the fact that $\mathbf{N}_{y}$ and $\mathbf{N}_{y}^{'}$ are iid. The expression $\var\left[\mathbf{N}_{y}f\left(\mathbf{N}_{y}\right)\right]$ is simply a sum of terms of the form of $\var\left[\mathbf{N}_{y}^{\gamma}\right]$ where $0<\gamma\leq1$. Even the covariance terms can be bounded by the square root of the product of these terms by the Cauchy Schwarz inequality.

Let $p_{y}=f_{Y}(y)$. Consider the Taylor series expansion of the function $h(x)=x^{\gamma}$ at the point $Np_{y}$. This is $$\begin{aligned}
h(x) & = & \left(Np_{y}\right)^{\gamma}+\gamma\left(Np_{y}\right)^{\gamma-1}\left(x-Np_{y}\right)+\frac{\gamma(\gamma-1)}{2}\left(Np_{y}\right)^{\gamma-2}\left(x-Np_{y}\right)^{2}\nonumber \\
 &  & +\sum_{k=3}^{\infty}\frac{\gamma(\gamma-1)\dots(\gamma-k+1)}{k!}\left(Np_{y}\right)^{\gamma-k}\left(x-Np_{y}\right)^{2}.\label{eq:taylor_gamma}\end{aligned}$$ From Riordan [@riordan1937moment], we know that the $i$th central moment of $\mathbf{N}_{y}$ is $O\left(N^{\left\lfloor i/2\right\rfloor }\right)$. Then since $\gamma\leq1$, the last terms in (\[eq:taylor\_gamma\]) are $O\left(N^{-1}\right)$ when $x=\mathbf{N}_{y}$ and we take the expectation. Thus $$\begin{aligned}
\bE\left[\mathbf{N}_{y}^{\gamma}\right] & = & \left(Np_{y}\right)^{\gamma}+\frac{\gamma(\gamma-1)}{2}\left(Np_{y}\right)^{\gamma-1}(1-p_{y})+O\left(N^{-1}\right)\\
\implies\bE\left[\mathbf{N}_{y}^{\gamma}\right]^{2} & = & \left(Np_{y}\right)^{2\gamma}+\gamma(\gamma-1)(1-p_{y})\left(Np_{y}\right)^{2\gamma-1}+\left(\frac{\gamma(\gamma-1)}{2}\right)^{2}\left(Np_{y}\right)^{2\gamma-2}\\
 &  & +O\left(N^{-1}\right).\end{aligned}$$ By a similar Taylor series expansion, we have that $$\bE\left[\mathbf{N}_{y}^{2\gamma}\right]=\left(Np_{y}\right)^{2\gamma}+\gamma(2\gamma-1)(1-p_{y})\left(Np_{y}\right)^{2\gamma-1}+O\left(N^{-1}\right).$$ Combining these results gives $$\begin{aligned}
\var\left[\mathbf{N}_{y}^{\gamma}\right] & = & \bE\left[\mathbf{N}_{y}^{2\gamma}\right]-\bE\left[\mathbf{N}_{y}^{\gamma}\right]^{2}\\
 & = & O\left(N^{2\gamma-1}+N^{2\gamma-2}+N^{-1}\right)\\
 & = & O\left(N\right),\end{aligned}$$ where the last step follows from the fact that $\gamma\leq1$. Combining this result with (\[eq:var\_cond\]) gives $$\var\left[\bE\left[\left.\g{h_{X},h_{X|Y}}\right|\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}\right]\right]=O\left(\frac{1}{N}\right).$$ By the law of total variance, $\var\left[\g{h_{X},h_{X|Y}}\right]=O\left(\frac{1}{N}\right)$.

Extension to the Generalized Case {#subsec:generalCase}
---------------------------------

In this section, we sketch the theory for the case where both $\X$ and $\Y$ have a mixture of discrete and continuous components. Denote the discrete and continuous components of $\mathbf{X}$ as $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, respectively. Similarly, denote the discrete and continuous components of $\Y$ as $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$, respectively. Then the generalized mutual information is $$\begin{aligned}
G(\X;\Y) & =\sum_{\substack{y_{1}\in\mathcal{S}_{Y_{1}}\\
x_{1}\in\mathcal{S}_{X_{1}}
}
}\int g\left(\frac{f_{X}(x_{1},x_{2})f_{Y}(y_{1},y_{2})}{f_{XY}(x_{1},x_{2},y_{1},y_{2})}\right)f_{XY}(x_{1},x_{2},y_{1},y_{2})dx_{2}dy_{2}\\
 & =\sum_{\substack{y_{1}\in\mathcal{S}_{Y_{1}}\\
x_{1}\in\mathcal{S}_{X_{1}}
}
}f_{X_{1}Y_{1}}(x_{1},y_{1})\int g\left(\frac{f_{X_{1}}(x_{1})f_{Y_{1}}(y_{1})f_{X_{2}|X_{1}}(x_{2}|x_{1})f_{Y_{2}|Y_{1}}(y_{2}|y_{1})}{f_{X_{1}Y_{1}}(x_{1},y_{1})f_{X_{2}Y_{2}|X_{1}Y_{1}}(x_{2},y_{2}|x_{1}y_{1})}\right)f_{X_{2}Y_{2}|X_{1}Y_{1}}(x_{2},y_{2}|x_{1}y_{1})dx_{2}dy_{2}.\end{aligned}$$ Define $\mathbf{N}_{y_{1}}=\sum_{i=1}^{N}1_{\{\mathbf{Y}_{1,i}=y_{1}\}}$ where $\mathbf{Y}_{1,i}$ is the discrete component of $\mathbf{Y}_{i}$. Then the estimator we use for $f_{Y_{1}}(y_{1})$ is $\mathbf{N}_{y_{1}}/N$. The estimators for $f_{X_{1}}(x_{1})$ and $f_{X_{1}Y_{1}}(x_{1},y_{1})$ are defined similarly with $\N_{x_{1}}$ and $\N_{z_{1}}$.

We first consider the conditional bias of the resulting plug-in estimator where we condition on the discrete random variables. Recall that by Taylor series expansions, we decompose the bias into “variance-like” terms in (\[eq:varterms\]) and “bias-like” terms in (\[eq:biasterms\]). For the bias-like term, if we condition on the discrete random variables, then the equivalent expression in (\[eq:E\_fx\]) in this case is multiplied by $\N_{x_{1}}/N$. This results in terms of the form of, for example, $\left(\frac{\N_{z_{1}}}{N}-f_{Z_{1}}(x_{1},y_{1})\right)^{i}$. The expected value of these expressions is the $i$th central moment of a binomial random variable divided by $N^{i}$ which is $O(1/N)$ for $i\geq1$. Thus these terms contribute $O(1/N)$ to the bias. In all other cases, the expected value of $\left(\frac{\N_{z_{1}}}{N}\right)^{i}$ is $O(1)$. Thus only the constants are affected by these terms in the equivalent expression in (\[eq:varterm\_bound\]). Similar results hold for the estimators of $f_{X_{1}}(x_{1})$ and $f_{Y_{1}}(y_{1})$.

For the “variance-like” terms, we can simply factor out the estimators for $f_{X_{1}}(x_{1})$, $f_{Y_{1}}(y_{1})$, and $f_{Z_{1}}(z_{1})$. The expected value of these estimators is again $O(1)$ so they only affect the constants.

For the variance, the law of total variance can again be used by conditioning on the discrete components. For the conditional variance, the Lipschitz conditions on $g$ in this case simply scales the resulting terms by the square of the estimators for $f_{X_{1}}(x_{1})$, $f_{Y_{1}}(y_{1})$, and $f_{Z_{1}}(z_{1})$. Then since the expected value of the square of these estimators is $O(1),$ the expected value of the conditional variance is still $O(1/N)$. Then by similar arguments given above for the bias and in Section \[subsec:varproof\], the variance of the conditional expectation of the estimator is also $O(1/N)$. Thus the total variance is $O(1/N)$.

Proof of Theorem 6 (CLT) {#sec:cltProof}
========================

This proof shares some similarities with the CLT proof for the divergence functional estimators in [@moon2016arxiv; @moon2016isit]. The primary differences again deal with handling products of marginal density estimators and with handling two of the terms in the Efron-Stein inequality. We will first find the asymptotic distribution of $$\begin{aligned}
\sqrt{N}\left(\gt-\bE\left[\gt\right]\right) & =\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-\ez{Z_{i}}\left[g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right)\\
 & +\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(\ez{Z_{i}}\left[g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]-\bE\left[g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right).\end{aligned}$$ By the standard central limit theorem [@durrett2010probability], the second term converges in distribution to a Gaussian random variable with variance $$\var\left[\ez Z\left[g\left(\frac{\ft X(\mathbf{X})\ft Y(\mathbf{Y})}{\ft Z(\mathbf{X},\mathbf{Y})}\right)\right]\right].$$ All that remains is to show that the first term converges in probability to zero as Slutsky’s theorem [@gut2012probability] can then be applied. Denote this first term as $\W_{N}$ and note that $\bE\left[\W_{N}\right]=0$.

We will use Chebyshev’s inequality combined with the Efron-Stein inequality to bound the variance of $\W_{N}$. Consider the samples $\left\{ \Z_{1},\dots,\Z_{N}\right\} $ and $\left\{ \Z_{1}^{'},\Z_{2},\dots,\Z_{N}\right\} $ and the respective sequences $\W_{N}$ and $\W_{N}^{'}$. This gives $$\begin{aligned}
\W_{N}-\W_{N}^{'} & =\frac{1}{\sqrt{N}}\left(g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-\ez{Z_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]\right)\nonumber \\
 & +\frac{1}{\sqrt{N}}\left(g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)-\ez{Z_{1}^{'}}\left[g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)\right]\right)\nonumber \\
 & +\frac{1}{\sqrt{N}}\sum_{i=2}^{N}\left(g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right).\label{eq:Wndiff}\end{aligned}$$

Note that $$\bE\left[\left(g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-\ez{Z_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]\right)^{2}\right]=\bE\left[\var_{\X_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]\right].$$ We will use the Efron-Stein inequality to bound $\var_{\X_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]$. We thus need to bound the conditional expectation of the term $$\left|g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{1})\ft Y^{'}(\mathbf{Y}_{1})}{\ft Z^{'}(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right|^{2},$$ where $\Z_{i}$ is replaced with $\Z_{i}^{'}$ in the KDEs for some $i\neq1$. Using similar steps as in Section \[sec:VarProof\], we have that $$\bE\left[\left|g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{1})\ft Y^{'}(\mathbf{Y}_{1})}{\ft Z^{'}(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right|^{2}\right]=O\left(\frac{1}{N^{2}}\right).$$ Then by the Efron-Stein inequality, $\var_{\X_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]=O\left(\frac{1}{N}\right)$. Therefore $$\bE\left[\frac{1}{N}\left(g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-\ez{Z_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]\right)^{2}\right]=O\left(\frac{1}{N^{2}}\right).$$ A similar result holds for the $g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)$ term in (\[eq:Wndiff\]).

For the third term in (\[eq:Wndiff\]), $$\begin{aligned}
 & \lefteqn{\bE\left[\left(\sum_{i=2}^{N}\left|g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right|\right)^{2}\right]}\\
 & =\sum_{i,j=2}^{N}\bE\left[\left|g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right|\right.\\
 & \times\left.\left|g\left(\frac{\ft X(\mathbf{X}_{j})\ft Y(\mathbf{Y}_{j})}{\ft Z(\mathbf{X}_{j},\mathbf{Y}_{j})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{j})\ft Y^{'}(\mathbf{Y}_{j})}{\ft Z^{'}(\mathbf{X}_{j},\mathbf{Y}_{j})}\right)\right|\right]\end{aligned}$$ For the $N-1$ terms where $i=j$, we know from Section \[sec:VarProof\] that $$\bE\left[\left|g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right|^{2}\right]=O\left(\frac{1}{N^{2}}\right).$$ Thus these terms contribute $O(1/N)$. For the $N^{2}-N$ terms where $i\neq j$, we can do multiple substitutions of the form $\mathbf{u}_{j}=\frac{\X_{j}-\X_{1}}{h_{X}}$ resulting in $$\begin{aligned}
\bE\left[\left|g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right|\right.\\
\times\left.\left|g\left(\frac{\ft X(\mathbf{X}_{j})\ft Y(\mathbf{Y}_{j})}{\ft Z(\mathbf{X}_{j},\mathbf{Y}_{j})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{j})\ft Y^{'}(\mathbf{Y}_{j})}{\ft Z^{'}(\mathbf{X}_{j},\mathbf{Y}_{j})}\right)\right|\right] & =O\left(\frac{h_{X}^{2d_{X}}h_{Y}^{2d_{Y}}}{N^{2}}\right).\end{aligned}$$ Since $h_{X}^{d_{X}}h_{Y}^{d_{Y}}=o(1)$, $$\bE\left[\left(\sum_{i=2}^{N}\left|g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right|\right)^{2}\right]=o(1).$$ Combining all of these results with Jensen’s inequality gives $$\begin{aligned}
\bE\left[\left(\W_{N}-\W_{N}^{'}\right)^{2}\right] & \leq\frac{3}{N}\bE\left[\left(g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)-\ez{Z_{1}}\left[g\left(\frac{\ft X(\mathbf{X}_{1})\ft Y(\mathbf{Y}_{1})}{\ft Z(\mathbf{X}_{1},\mathbf{Y}_{1})}\right)\right]\right)^{2}\right]\\
 & +\frac{3}{N}\bE\left[\left(g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)-\ez{Z_{1}^{'}}\left[g\left(\frac{\ft X(\mathbf{X}_{1}^{'})\ft Y(\mathbf{Y}_{1}^{'})}{\ft Z(\mathbf{X}_{1}^{'},\mathbf{Y}_{1}^{'})}\right)\right]\right)^{2}\right]\\
 & +\frac{3}{N}\bE\left[\left(\sum_{i=2}^{N}\left(g\left(\frac{\ft X(\mathbf{X}_{i})\ft Y(\mathbf{Y}_{i})}{\ft Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)-g\left(\frac{\ft X^{'}(\mathbf{X}_{i})\ft Y^{'}(\mathbf{Y}_{i})}{\ft Z^{'}(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right)\right)^{2}\right]\\
 & =o\left(\frac{1}{N}\right).\end{aligned}$$ Applying the Efron-Stein inequality gives that $\var\left[\W_{N}\right]=o(1)$. Then by ChebyShev’s inequality, $\W_{N}$ converges to zero in probability. This completes the proof for the plug-in estimator.

For the weighted ensemble estimator, we can write $$\begin{aligned}
\sqrt{N}\left(\g w-\bE\left[\g w\right]\right) & =\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})\left(g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right.\\
 & \left.-\ez{Z_{i}}\left[g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right)\\
 & +\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(\ez{Z_{i}}\left[\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right.\\
 & \left.-\bE\left[\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right).\end{aligned}$$ By the central limit theorem, the second term converges in distribution to a zero-mean Gaussian random variable with variance $$\var\left[\ez{Z_{i}}\left[\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right].$$ From the previous results, the first term converges to zero in probability as it can be written as $$\begin{aligned}
\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\left(g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right.\\
\left.-\ez{Z_{i}}\left[g\left(\frac{\ftl X(\mathbf{X}_{i})\ftl Y(\mathbf{Y}_{i})}{\ftl Z(\mathbf{X}_{i},\mathbf{Y}_{i})}\right)\right]\right) & =\sum_{l_{X}\in\mathcal{L}_{X},l_{Y}\in\mathcal{L}_{Y}}w(l_{X},l_{Y})o_{P}(1)\\
 & =o_{P}(1),\end{aligned}$$ 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 7 follows a similar procedure as the extension to the ensemble case.

[^1]: A.A@university.edu

[^2]: B.B@university.edu

[^3]: D.D@university.edu
