Ensemble Estimation of Generalized Mutual Information with Applications to Genomics

TL;DR

This paper introduces GENIE, an ensemble estimator for generalized mutual information that achieves parametric convergence rates in mixed discrete-continuous data, with applications in genomics and other fields.

Contribution

The paper develops the first nonparametric mutual information estimator that attains parametric rates for mixed data, using an ensemble approach with simple implementation.

Findings

Achieves $1/N$ mean squared error convergence rate in continuous and mixed cases.

Demonstrates effectiveness on simulated data and gene relationship analysis in genomics.

Provides theoretical guarantees including a central limit theorem and minimax rates.

Abstract

Mutual information is a measure of the dependence between random variables that has been used successfully in myriad applications in many fields. Generalized mutual information measures that go beyond classical Shannon mutual information have also received much interest in these applications. We derive the mean squared error convergence rates of kernel density-based plug-in estimators of general mutual information measures between two multidimensional random variables $\mathbf{X}$ and $\mathbf{Y}$ for two cases: 1) $\mathbf{X}$ and $\mathbf{Y}$ are continuous; 2) $\mathbf{X}$ and $\mathbf{Y}$ may have any mixture of discrete and continuous components. Using the derived rates, we propose an ensemble estimator of these information measures called GENIE by taking a weighted sum of the plug-in estimators with varied bandwidths. The resulting ensemble estimators achieve the $1/N$ parametric mean squared error convergence rate when the conditional densities of the continuous variables are sufficiently smooth. To the best of our knowledge, this is the first nonparametric mutual information estimator known to achieve the parametric convergence rate for the mixture case, which frequently arises in applications (e.g. variable selection in classification). The estimator is simple to implement and it uses the solution to an offline convex optimization problem and simple plug-in estimators. A central limit theorem is also derived for the ensemble estimators and minimax rates are derived for the continuous case. We demonstrate the ensemble estimator for the mixed case on simulated data and apply the proposed estimator to analyze gene relationships in single cell data.

Tables

Table 1. Table I: Negative log-log slope of the MSE for the plots in Figure  1 as a function of the sample size. The estimated values are compared to the theoretical values. While the estimators outperform the theoretical rates in this finite sample regime, the estimator rates generally approach the theoretical rates as the dimension increases.

Estimator	$\left|S_Y\right|=3,d=4$	$\left|S_Y\right|=3,d=6$	$\left|S_Y\right|=3,d=9$	$\left|S_Y\right|=6,d=6$
Kernel Estimator	0.90	0.59	0.32	0.49
Kernel Theoretical	0.50	0.33	0.22	0.33
Weighted Estimator	1.45	1.42	1.21	1.51
Weighted Theoretical	1	1	1	1

Table 2. Table II: Estimated Rényi MI and DREMI between collections of genes and cell type for mouse bone marrow scRNA-seq data  [ 94 ] and the corresponding classification accuracies from a linear support vector machine and random forests using 10-fold cross validation. Gene collections are selected from the KEGG pathways associated with the hematopoietic cell lineage. The fifth column (with heading “Combined”) gives the result when combining all genes together from the four KEGG pathways. The last column gives the results when selecting 8 genes at random averaged over 50 trials. MI results are presented in the form of mean ± plus-or-minus \pm standard deviation which are calculated from 1000 bootstrapped samples.

	Platelets	Erythrocytes	Neutrophils	Macrophages	Combined	Random
Estimated Rényi MI	$0.24\pm 0.11$	$0.66\pm 0.11$	$0.27\pm 0.10$	$0.15\pm 0.09$	$1.65\pm 0.36$	$0.007\pm 0.07$
Estimated DREMI	$0.25\pm 0.22$	$0.04\pm 0.03$	$0.20\pm 0.12$	$0.41\pm 0.45$	$0.88\pm 0.35$	$0.03\pm 0.08$
SVM Accuracy	57.4%	57.5%	52.9%	52.9%	65.4%	43.2%
Random Forests Accuracy	60.3%	60.0%	57.8%	57.8%	65.9%	52.3%

Table 3. Table III: Results when computing the Rényi MI, DREMI, and the SIS between the time course variable and a single gene (i.e. d X = 1 subscript 𝑑 𝑋 1 d_{X}=1 ) for all genes in the data with nonnegative correlations with time. The top 10 genes for each branch and score are shown here. The SIS score corresponds to the correlation coefficient. Many of the genes are known to be associated with their respective tissues.

Neural Progenitors Branch						Neural Crest Branch
Rényi MI		DREMI		SIS		Rényi MI		DREMI		SIS
LINC00526	1.004	BRWD1-AS2	8.768	FOS	0.934	HOXB7	0.837	CRYL1	10.324	RARB	0.918
GTF2E2	1.003	NR2E1	8.579	GTF2E2	0.934	SEPT6	0.820	SHC4	9.890	DDIT4	0.917
SEPT6	0.966	ZNF804A	8.505	SLC18B1	0.932	HOXA3	0.818	SLITRK2	9.692	HOXB7	0.909
SREBF2	0.963	SYT4	8.233	JAM2	0.931	HOXA7	0.818	GDNF-AS1	9.304	RGCC	0.908
EFCAB1	0.948	NTNG1	8.146	EGR1	0.930	RBP1	0.806	PRAME	9.235	IGFBP7	0.905
RP11-68606.2	0.937	GPR1	8.001	CX3CL1	0.929	ACADS	0.804	PAQR6	9.164	HOXA5	0.904
B2M	0.936	POU3F4	7.899	MAGEL2	0.927	HOXB5	0.803	Clorf198	9.044	AEBP1	0.903
CX3CL1	0.928	HSD17B8	7.704	LINC00632	0.927	HOXB6	0.794	HSPB2	8.971	HOXA7	0.901
RP3-525N10.2	0.928	LINC00092	7.686	TPPP3	0.927	HOXB3	0.793	LINC00518	8.904	ACADS	0.901
C20orf96	0.926	WNT4	7.685	GSTM3	0.927	RND3	0.791	AZGP1	8.890	PPP1R15A	0.900

Table 4. Table IV: Results when computing the Rényi MI, DREMI, and the SIS between the time course variable and multiple genes using a greedy forward-selection approach. The top 10 genes for each branch and score are shown here. The SIS score corresponds to the correlation coefficient of gene expression with the regression residuals. Many of these genes are known to be associated with their respective tissues.

Neural Progenitors (NP) Branch						Neural Crest (NC) Branch
Rényi MI		DREMI		SIS		Rényi MI		DREMI		SIS
LINC00526	1.004	BRWD1-AS2	9	FOS	0.934	HOXB7	0.837	CRYL1	10	RARB	0.918
HFE	1.382	FOSL1	15	CDC37L1-AS1	0.360	AF127936.9	1.119	IDH3B	14	ID3	0.394
DOC2A	1.675	TCP11	30	SH3GL2	0.194	RGR	1.407	AC142528.1	11	ZNF564	0.282
P4HA1	1.931	BRD9	48	ATP13A3	0.229	RP11-324E6.10	2.151	CFL1	33	MALRD1	0.188
HIST1H1C	2.030	HFE	103	DARS	0.178	LIMA1	2.347	RP11-676J15.1	64	BRWD1-AS2	0.165
PRDM12	2.152	RP11-225H22.4	224	RP11-390P2.4	0.171	ITGA9	2.581	RP5-1098D14.1	87	RP11-10A14.4	0.169
ACTR3C	2.189	ATG9A	335	ZNF484	0.167	TTLL9	2.699	BMP8B	193	SLC10A5	0.126
TRDC	2.222	RP11-35015.1	848	SIDT2	0.164	ITPKB	2.826	SLITRK2	936	RP3-402G11.26	0.117
MMEL1	2.285	GAS2L3	2385	NEFH	0.140	ACOT1	2.857	TMCC	1846	ABCC1	0.122
LINC01229	2.327	PRAC1	5845	ZNF587	0.127	DPYSL4	2.869	MLANA	6201	RPSA	0.124

Table 5. Table V: Results when computing the Rényi MI, DREMI, and the SIS between the time course variable and multiple genes using a greedy forward-selection approach when starting with three or four relevant genes identified in  [ 38 ] . The top 10 genes for each branch and score are shown here. The SIS score corresponds to the correlation coefficient of gene expression with the regression residuals. Many of these genes are known to be associated with their respective tissues.

Neural Progenitors (NP) Branch						Neural Crest (NC) Branch
Rényi MI		DREMI		SIS		Rényi MI		DREMI		SIS
NKX2-8	-	NKX2-8	-	NKX2-8	-	PAX3	-	PAX3	-	PAX3	-
EN2	-	EN2	-	EN2	-	FOXD3	-	FOXD3	-	FOXD3	-
SOX1	-	SOX1	-	SOX1	-	SOX9	-	SOX9	-	SOX9	-
SLX1B	2.030	MSS51	15	FAM174A	0.567	SOX10	-	SOX10	-	SOX10	-
PTK6	2.613	KNSTRN	247	RAB27A	0.329	RP11-867G23.8	1.614	GRHPR	7	ARMC5	0.543
SIRT3	3.204	TEKT3	548	CTD-2568A17.1	0.321	ZBED3	1.759	LINC01389	23	ABL1	0.325
HEY1	3.596	FAM72C	1998	CECR5	0.209	RP3-460G2.2	1.830	CTB-25B13.5	50	PSPN	0.338
OLIG2	4.021	GZMK	4904	LOXL1	0.163	LRSAM1	1.855	CFAP58	68	RP11-354P11.3	0.262
MZF1-AS1	4.689	QRICH2	16337	PROB1	0.147	NOTUM	1.930	ZNF774	271	WWTR1-AS1	0.185
BMPR1B	5.133	PPFIA4	42236	YY2	0.161	CPLX2	1.945	LINC00327	588	KIF25	0.144

Equations

$$D_g(P||Q):={\mathbb{E}}_Q\left[g\left(\frac{dP}{dQ}\right)\right],$$

$$I(\mathbf{X};\mathbf{Y}):=D_g\left(P_X{P}_Y\Vert P_{XY}\right).$$

$$I(\mathbf{X};\mathbf{Y})=\int g\left(\frac{f_X\left(x\right)f_Y\left(y\right)}{f_{XY}\left(x,y\right)}\right)f_{XY}\left(x,y\right)dxdy.$$

$$I(\mathbf{X};\mathbf{Y})$$

$$=\sum _{z_D}{f}_{Z_D}\left(z_D\right)\int g\left(R_1\left(z_C\right)R_2\left(z_D\right)\right)f_{Z_C|Z_D}\left(z_C|z_D\right)d{z}_C,$$

$$R_1(z_C)$$

$$R_2\left(z_D\right)$$

$$I(\mathbf{X};\mathbf{Y})=\int g\left(\frac{f_X\left(x\right)f_Y\left(y\right)}{f_{XY}\left(x,y\right)}\right)f_{XY}\left(x,y\right)dxdy.$$

$$I_{\nu }(\mathbf{X};\mathbf{Y})=\int g\left(\frac{f_X\left(x\right)f_Y\left(y\right){\nu }_1{\nu }_2}{f_{XY}\left(x,y\right){\nu }_3}\right)f_{XY}\left(x,y\right)dxdy.$$

$${\tilde{\mathbf{f}}}_{X,h_X}(x)$$

$${\tilde{\mathbf{f}}}_{Y,h_Y}(y)$$

$${\tilde{\mathbf{f}}}_{Z,h_Z}(x,y)$$

$${\tilde{\mathbf{G}}}_{h_X,h_Y}=\frac{1}{N}\sum _{i=1}^{N}g\left(\frac{{\tilde{\mathbf{f}}}_{X,h_X}({\mathbf{X}}_i){\tilde{\mathbf{f}}}_{Y,h_Y}({\mathbf{Y}}_i){\nu }_1{\nu }_2}{{\tilde{\mathbf{f}}}_{Z,h_Z}({\mathbf{X}}_i,{\mathbf{Y}}_i){\nu }_3}\right).$$

$$\left|D^{q}f(x)-D^{q}f(y)\right|\le H{\Vert x-y\Vert }^{s-|q|},$$

$${\int }_{x\in \mathcal{S}}{\left({\int }_{u:K(u)>0,x+uh\notin \mathcal{S}}{K}^l(u)u^q{D}^{q}f(x)du\right)}^{t}dx=v_t(h),$$

$$v_t(h)=\sum _{i=1}^{r-|q|}{e}_{i,q,t,l}{h}^i+o\left(h^{r-|q|}\right),$$

$$\mathbb{B}\left[{\tilde{\mathbf{G}}}_{h_X,h_Y}\right]$$

$$\mathbb{V}\left[{\tilde{\mathbf{G}}}_{h_X,h_Y}\right]\le \frac{22C_g^2||K_X\cdot K_Y|{|}_{\infty }^2}{N}.$$

$$\underset{N\to \infty }{\mathrm{liminf}}\underset{{\hat{G}}_N}{\inf }\underset{f_{XY}\in \Sigma \left(s,H,{ϵ}_0,{ϵ}_{\infty }\right)}{\sup }\mathbb{E}\left[{\left({\hat{G}}_N-I(\mathbf{X};\mathbf{Y})\right)}^2\right]\ge c{N}^{-\gamma }.$$

$${\mathbf{N}}_y$$

$${\mathbf{N}}_x$$

$${\mathbf{N}}_{xy}$$

$${\mathcal{X}}_x$$

$${\mathcal{Y}}_y$$

$${\tilde{\mathbf{f}}}_{X_C|x,h_{X_C|x}}\left(x\right)$$

$${\tilde{\mathbf{f}}}_{Y_C|y,h_{Y_C|y}}\left(y\right)$$

$${\tilde{\mathbf{f}}}_{Z_C|z,h_{Z_C|z}}\left(x,y\right)$$

$${\tilde{\mathbf{G}}}_{h_{X_C|x},h_{Y_C|y}}=\frac{1}{{\mathbf{N}}_{xy}}\sum _{{\mathbf{X}}_C\in {\mathcal{X}}_x\text{AND}{\mathbf{Y}}_C\in {\mathcal{Y}}_y}g\left(\frac{{\tilde{\mathbf{f}}}_{X_C|x,h_{X_C|x}}\left({\mathbf{X}}_C\right){\tilde{\mathbf{f}}}_{Y_C|y,h_{Y_C|y}}\left({\mathbf{Y}}_C\right)}{{\tilde{\mathbf{f}}}_{Z_C|z,h_{Z_C|z}}\left({\mathbf{X}}_C,{\mathbf{Y}}_C\right)}\times \frac{{\mathbf{N}}_x{\mathbf{N}}_y}{N{\mathbf{N}}_{xy}}\right).$$

$${\tilde{\mathbf{G}}}_{h_{X_C|X_D},h_{Y_C|Y_D}}=\sum _{x\in {\mathcal{S}}_{X_D},y\in {\mathcal{S}}_{Y_D}}\frac{{\mathbf{N}}_{xy}}{N}{\tilde{\mathbf{G}}}_{h_{X_C|x},h_{Y_C|y}}.$$

$$\mathbb{E}\left[{\mathbf{N}}_{xy}^{\alpha }\right]$$

$$\mathbb{E}\left[{\mathbf{N}}_{xy}^{\lambda }{\mathbf{N}}_x^{\beta }{\mathbf{N}}_y^{\gamma }\right]$$

$$\mathbb{B}\left[{\tilde{\mathbf{G}}}_{h_{X_C|X_D},h_{Y_C|Y_D}}\right]$$

$$s\alpha$$

$$s\beta$$

$$1-\beta d_X-\alpha d_Y$$

$${\tilde{\mathbf{G}}}_w=\sum _{l\in \mathcal{L}}w(l)\sum _{x\in {\mathcal{S}}_{X_D},y\in {\mathcal{S}}_{Y_D}}\frac{{\mathbf{N}}_{xy}}{N}{\tilde{\mathbf{G}}}_{h_{X_C|x}(l),h_{Y_C|y}(l)}.$$

Full text

1. $\mathbf{X}$ and $\mathbf{Y}$ are continuous; 2) $\mathbf{X}$ and $\mathbf{Y}$ may have a mixture of discrete and continuous components. Using the derived rates, we propose an ensemble estimator of these information measures called GENIE by taking a weighted sum of the plug-in estimators with varied bandwidths. The resulting ensemble estimators achieve the $1/N$ parametric mean squared error convergence rate when the conditional densities of the continuous variables are sufficiently smooth. To the best of our knowledge, this is the first nonparametric mutual information estimator known to achieve the parametric convergence rate for the mixture case, which frequently arises in applications (e.g. variable selection in classification). The estimator is simple to implement and it uses the solution to an offline convex optimization problem and simple plug-in estimators. A central limit theorem is also derived for the ensemble estimators and minimax rates are derived for the continuous case. We demonstrate the ensemble estimator for the mixed case on simulated data and apply the proposed estimator to analyze gene relationships in single cell data.

I Introduction

Mutual information (MI) is a measure of the amount of shared information between a pair of random variables $\mathbf{X}$ and $\mathbf{Y}$. MI estimation is related to the problem of estimating functionals of probability distributions, which has received deserved attention in recent years [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Many statistical problems rely in some form upon accurate estimation of functionals of probability distributions including estimating the decay rates of error probabilities [18], estimating bounds on the Bayes error rate [19, 20, 21, 22], and hypothesis testing [23, 24, 25]. MI estimation, in particular, also has many applications in information theory and machine learning including independent subspace analysis [26], structure learning [27], fMRI data processing [28], forest density estimation [29], clustering [30], neuron classification [31], blind source separation [32], intrinsically motivated reinforcement learning [33, 34], as well as other data science applications such as sociology [35], computational biology [36, 37, 38], and improving neural network models [39]. A particularly common application is feature selection or extraction where features are chosen to maximize the MI between the chosen features (represented by $\mathbf{X}$) and the outcome variables (represented by $\mathbf{Y}$) [40, 41, 42, 43, 44].

In many of these applications, the variables $\mathbf{X}$ and $\mathbf{Y}$ may have any mixture of discrete and continuous components. In feature selection, for example, the predictor labels may have discrete components (e.g. classification labels) while the input variables may have a mixture of discrete and continuous features. To the best of our knowledge, there are currently no nonparametric MI estimators that are known to achieve the parametric mean squared error (MSE) convergence rate $1/N$ ($N$ is the number of samples) in this setting where $\mathbf{X}$ and/or $\mathbf{Y}$ contain a mixture of discrete and continuous components. Instead, most existing estimators of MI focus on the cases where both $\mathbf{X}$ and $\mathbf{Y}$ are either purely discrete or purely continuous. Also, while many nonparametric estimators of MI exist, most have not been generalized beyond Shannon or Rényi information. Furthermore, minimax convergence rates are currently unknown for the continuous and the mixture cases.

In this paper, we provide a framework for nonparametric estimation of a large class of MI measures where we only have available a finite population of i.i.d. samples.** **This framework can be applied to accurately estimate general MI measures when either $\mathbf{X}$ and $\mathbf{Y}$ are purely continuous or the mixed case when $\mathbf{X}$ and $\mathbf{Y}$ may contain a mixture of discrete and continuous components. We derive an MI estimator for these cases that achieves the parametric MSE rate when the conditional densities of the continuous variables are sufficiently smooth, thus achieving the minimax rate (which we also derive) in this setting. We call this estimator the Generalized ENsemble Information Estimator (GENIE).

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

I-A Related Work

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

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

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

Recent work has focused on nonparametric divergence estimation for continuous random variables. One approach [7, 8, 9, 10] uses an optimal KDE to achieve the parametric convergence rate when the densities are at least $d$ [9, 10] or $d/2$ [7, 8] times differentiable where $d$ is the dimension of the data. These methods, like ours, assume that the densities are bounded away from zero as this simplifies the analysis. However, this induces a boundary on the densities’ support set. For accurate estimation, the optimal KDE approaches require knowledge of the density support boundary and are difficult to construct near the boundary. Numerical integration may also be required for estimating some divergence functionals under this approach, which can be computationally expensive. In contrast, our approach to MI estimation does not require numerical integration and can be performed without knowledge of the support boundary.

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

More closely related work [11, 13, 12, 17, 14, 64, 15, 16] uses an ensemble approach to estimate entropy or divergence functionals for continuous random variables. These works construct an ensemble of simple plug-in estimators by varying the neighborhood size of density estimators. They then take a weighted average of the estimators where the weights are chosen to decrease the bias with only a small increase in the variance. The parametric rate of convergence is achieved when the densities are either $d$ [11, 13, 12, 15] or $d/2$ [14, 64, 16] times differentiable. These approaches are simple to implement as they only require simple plug-in estimates and the solution of an offline convex optimization problem. The ensemble approach also automatically corrects for bias at the boundary of the densities’ support set.

Finally, [65] showed that $k$-nn or KDE based approaches underestimate the MI when the MI is large. As MI increases, the dependencies between random variables increase which results in less smooth densities. Thus this isn’t an issue when the densities are smooth [7, 8, 9, 10, 11, 13, 12, 14].

For the mixture setting, we focus on the important special case where the components of each observation are assumed to decompose into discrete and continuous dimensions. This enables the density to be factored: $f_{X}(x)=f_{X_{D}}\left(x_{D}\right)f_{X_{C}|X_{D}}\left(x_{C}|x_{D}\right)$ where $x_{C}$ and $x_{D}$ are the continuous and discrete components of $x$. We note that this excludes the more general case considered by [2] where one or more components can have discrete and continuous values simultaneously. However, our setting is a common occurrence in many machine learning and statistical problems. For example, a search within the UCI Machine Learning Repository [66] yields many datasets with such structure. Many statistical models have also focused on similar settings [67, 68, 69, 70]. Thus we believe that this special case warrants its own treatment and retain the more general case for future work. Despite the importance of this mixed setting, no other MI estimators have been derived or analyzed that achieve the parametric MSE convergence rate.

I-B Contributions

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

For purely continuous random variables, we derive the asymptotic bias and variance of kernel density plug-in MI estimators for general MI measures without boundary correction [71] (Section III). 2. 2.

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

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

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

We derive a central limit theorem for the ensemble estimators (Section V-B). 6. 6.

We apply the method to single-cell RNA-sequencing feature selection problems (Section VI).

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

II Mutual Information Functionals

We first define a family of MI functionals based on $f$-divergence functionals which are defined as follows. Let $P$ and $Q$ be probability measures on the Euclidean space $\mathcal{S}$. Let $g:(0,\infty)\rightarrow\mathbb{R}$ be the $f$-divergence shaping function. The $f$-divergence functional associated with $g$ is [77, 78]

[TABLE]

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

A generalized MI functional can be derived from (1). Let $\mathbf{X}$ and $\mathbf{Y}$ be (potentially multivariate) random variables with respective marginal probability measures $P_{X}$ and $P_{Y}$ and joint probability measure $P_{XY}$. Let $g$ be as before. Then the MI functional associated with $g$ is

[TABLE]

Shannon MI can be obtained from (2) by setting $g(t)=-\log t$.

If $\mathbf{X}$ and $\mathbf{Y}$ are purely continuous random variables with respective marginal probability densities $f_{X}$ and $f_{Y}$ and joint probability density $f_{XY}$, then (2) can be written as

[TABLE]

However, we are also interested in the case where $\mathbf{X}$ or $\mathbf{Y}$ may have a mixture of discrete and continuous components. In this special case, the distributions can be factored into a product of the conditional density and the probability mass functions. The MI can then be expressed as a sum of integrals which can then be individually estimated. To do this, denote the continuous and discrete components of $\mathbf{X}$ as $\mathbf{X}_{C}$ and $\mathbf{X}_{D}$, respectively. Denote $\mathbf{Y}_{C}$ and $\mathbf{Y}_{D}$ similarly. Let $\mathbf{Z}=\left(\mathbf{X},\mathbf{Y}\right)^{T}$ and let $\mathbf{Z}_{C}$ and $\mathbf{Z}_{D}$ be the respective continuous and discrete components of $\mathbf{Z}$. Consider the probability distributions $f_{XY}$, $f_{X}$, $f_{Y}$ and the corresponding densities that are obtained by conditioning on $\mathbf{X}_{D}$ and $\mathbf{Y}_{D}$, e.g. $f_{XY}\left(x_{C},x_{D},y_{C},y_{D}\right)=f_{X_{C}Y_{C}|X_{D}Y_{D}}\left(x_{C},y_{C}|x_{D},y_{D}\right)f_{X_{D}Y_{D}}\left(x_{D},y_{D}\right)$. Then after factoring the distributions, (2) can be written as

[TABLE]

where

[TABLE]

The expression $R_{1}$ is the ratio of the product of the conditional densities $f_{X_{C}|X_{D}}$ and $f_{Y_{C}|Y_{D}}$ to the conditional density $f_{X_{C}Y_{C}|X_{D}Y_{D}}$. It is a continuous function of $z_{C}$. Similarly, the expression $R_{2}$ is the ratio of the product of the probability mass functions (pmf) $f_{X_{D}}$ and $f_{Y_{D}}$ to the pmf $f_{X_{D}Y_{D}}$ and is a discrete function of $z_{D}$.

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

III Continuous Random Variables

For this section, we define KDE plug-in estimators of general MI measures under the assumption that $\mathbf{X}$ and $\mathbf{Y}$ are purely continuous. Thus $\mathbf{X}_{C}=\mathbf{X}$ and $\mathbf{Y}_{C}=\mathbf{Y}$ and we can write

[TABLE]

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

To more easily generalize our results to the mixture case, we consider a modified version of (5) where the densities are weighted as follows. Let $\nu$ be a 3-dimensional vector with $0<\nu_{i}\leq 1$ for each $i\in\{1,2,3\}$. We can then write

[TABLE]

The expression in (6) reduces to that in (5) when $\nu_{i}=1$ for each $i\in\{1,2,3\}$. When we generalize to the mixture case, the pmf estimators will be substituted into $\nu$.

III-A The KDE Plug-in Estimator

Let $f_{X}(x)$, $f_{Y}(y)$, and $f_{XY}(x,y)$ be $d_{X}$, $d_{Y}$, and $d_{X}+d_{Y}=d$-dimensional probability densities. Since we are assuming for now that $\mathbf{X}$ and $\mathbf{Y}$ are continuous with marginal densities $f_{X}$ and $f_{Y}$, the MI functional $I_{v}(\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 $h_{X}$, $h_{Y}$ be kernel bandwidths. Let $K_{X}(\cdot)$ and $K_{Y}(\cdot)$ be symmetric kernel functions with $\int K_{X}(x)dx=\int K_{Y}(y)dy=1$, $||K_{X}||_{\infty},\,||K_{Y}||_{\infty}<\infty$ where $||K||_{\infty}=\sup_{x}|K(x)|$. The KDEs for $f_{X}$, $f_{Y}$, and $f_{XY}=f_{Z}$, respectively, are

[TABLE]

where $h_{Z}=(h_{X},h_{Y})$. Then $I_{\nu}(\mathbf{X};\mathbf{Y})$ can be estimated with a KDE plug-in estimator:

[TABLE]

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

III-B MSE Convergence Rate of the Continuous Plug-in Estimator

We are interested in the MSE convergence rate of the KDE plug-in estimator in eq. (10). The MSE of an estimator can be expressed as the sum of the squared bias and the variance of the estimator. We first focus on the bias of the estimator $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$. The bias of nonparametric estimators typically depends on the smoothness of the functions that are being estimated. In our case, we have multiple functions including the joint and marginal densities and the function $g$. We quantify the smoothness of the densities using the Hölder class $\Sigma(s,H)$:

Definition 1** (Hölder Class).**

Let $\mathcal{X}\subset\mathbb{R}^{d}$ be a compact space. For $q=(q_{1},\dots,q_{d}),$ $q_{i}\in\mathbb{N},$ define $|q|=\sum_{i=1}^{d}q_{i}$ and $D^{q}=\frac{\partial^{|q|}}{\partial x_{1}^{q_{1}}\dots\partial x_{d}^{q_{d}}}$. The Hölder class $\Sigma(s,H)$ of functions on $L_{2}(\mathcal{X})$ consists of the functions $f$ that satisfy

[TABLE]

for all $x,\,y\in\mathcal{X}$ and for all $q$ s.t. $|q|\leq\left\lfloor s\right\rfloor$.

A key fact that comes from Definition 1 is that if a function $f$ belongs to $\Sigma(s,H),$ then it is $r=\left\lfloor s\right\rfloor$ times differentiable. Given this definition, the full assumptions we make to derive bias convergence rates 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\geq 2$.

• •

$(\mathcal{A}.3)$: $g(t_{1}/t_{2})$ 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,\dots$ 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}$, $f$ either $f_{X}$ or $f_{Y}$, and $d$ either $d_{X}$ or $d_{Y}$. Let $q=(q_{1},\dots,q_{d})$ with $q_{i}\in\mathbb{N}$ and $|q|\leq r=\left\lfloor s\right\rfloor$. Then we assume for any positive integers $t$ and $l$ that

[TABLE]

where $v_{t}(h)$ admits the expansion

[TABLE]

for some constants $e_{i,q,t,l}$.

These assumptions can largely be summarized as follows: 1) $f_{X}$, $f_{Y}$, $f_{XY}$, and $g$ are smooth ($\mathcal{A}.2$-$\mathcal{A}.4$) ; 2) $f_{X}$ and $f_{Y}$ have bounded support sets $\mathcal{S}_{X}$ and $\mathcal{S}_{Y}$ with respective dimensions $d_{X}$ and $d_{Y}$ ($\mathcal{A}.1$); 3) $f_{X}$, $f_{Y}$, and $f_{XY}$ are strictly lower bounded on their support sets ($\mathcal{A}.1$); and 4) the boundary of the support set is smooth ($\mathcal{A}.5$). More specifically, assumption $\mathcal{A}.5$ states that the support of the density is smooth with respect to the kernel $K$ in the sense that the expected value of a polynomial with coefficients consisting of the densities and their derivatives near the boundary is a smooth function of the bandwidth $h$. The inner integral in (11) captures this expectation while the outer integral averages this inner integral over all points near the boundary of the support. The $v_{t}(h)$ term captures the fact that the smoothness of this expectation is proportional to the smoothness of the function $D^{q}f(x)$.

While these assumptions may appear highly technical, they are satisfied for relatively simple support sets and for common kernels, functions $g$, and densities and thus are widely applicable [14, 64]. These assumptions are also comparable to those in similar studies on asymptotic convergence analysis [73, 12, 13, 10, 9, 11, 7, 8, 14]. Some studies consider the case where the densities are not strictly lower bounded, which makes the problem different [62, 63] with different minimax rates (see [63] for the entropy estimation case).

In particular, assumption $\mathcal{A}.5$ is satisfied if the kernel $K$ is smooth, has either circular or rectangular support (which includes product kernels), and the density support set consists of the unit cube. See Appendix C for details. The unit cube assumption is common in the nonparametric density functional estimation literature [56, 7, 9, 10, 11, 1] as the results can then be extended to density support sets that are isomorphic to the unit cube.

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

Proposition 1**.**

Let the density support set $\mathcal{S}$ be the unit cube. Let the derivatives of the density $f$ up to order $r$ vanish at the boundary of the density support set. Assume that $||K||_{\infty}<\infty$ and the support of $K$ is bounded with either rectangular or circular support. Then assumption $\mathcal{A}.5$ is satisfied.

The proof is given in Appendix C-C. The assumption of vanishing density derivatives at the boundary is strictly weaker than assumption $\mathcal{A}.5$. As an example, consider a standard Gaussian distribution truncated to the support $[-1,1]^{d}$. Clearly, the derivatives of this density do not vanish at the boundary. However, we show in Appendix C-D that this density satisfies $\mathcal{A}.5$.

We note that the boundary assumption $\mathcal{A}.5$ does not directly result in parametric convergence rates for the plug-in estimator $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$, which is in contrast with the boundary assumptions in [9, 10, 7, 8]. The estimators in [9, 10, 7, 8] 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 $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ without boundary correction, greatly simplifying our estimator.

Theorem 2** (Bias Expansion for Continuous $\mathbf{X},\mathbf{Y}$).**

Under assumptions $\mathcal{A}.0$-$\mathcal{A}.5$, the bias of $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ is

[TABLE]

where the constants in (12) are independent of the bandwidths $h_{X}$ and $h_{Y}$ and depend on the densities and their derivatives, the functional $g$ and its derivatives, and the kernels. They also include polynomial terms of $\nu_{1}\nu_{2}$ and $\nu_{3}$ when $\nu_{i}\neq 1$.

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

We also derive a refined expression for the bias that enables us to achieve the parametric convergence rate under less restrictive smoothness assumptions on the densities ($s>(d_{X}+d_{Y})/2$ compared to $s\geq d_{X}+d_{Y}$ for (12)). However, the resulting expansion has more terms and the ensemble estimator is more complicated to implement. Thus we have chosen to present the simpler case here. The more complex expansion and estimator are presented in Appendix B-B.

Having obtained an expression for the bias of $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$, we now present an upper bound on its variance to complete the derivation of its MSE.

Theorem 3** (Variance Bound for Continuous $\mathbf{X},\mathbf{Y}$).**

If the functional $g$ is Lipschitz continuous in both of its arguments with Lipschitz constant $C_{g}$, then the variance of $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ is

[TABLE]

The Lipschitz assumption on $g$ for the variance result is comparable to assumptions made by others for nonparametric estimation of distributional functionals [8, 9, 10, 72, 7] and is satisfied for Shannon and Renyi informations when the densities are bounded above and below. Note that Theorem 3 requires much less strict assumptions than Theorem 2. The proofs of Theorems 2 and 3 are given in Appendix D and E, respectively.

Theorems 2 and 3 indicate that for the MSE to go to zero, we require $h_{X},h_{Y}\rightarrow 0$ and $Nh_{X}^{d_{X}}h_{Y}^{d_{Y}}\rightarrow\infty$. In Section IV, we will use Theorems 2 and 3 to derive bias and variance expressions for the MI plug-in estimators under the more general cases where $\mathbf{X}$ and/or $\mathbf{Y}$ may contain a mixture of discrete and continuous components. We will then use these convergence rate results to derive MI ensemble estimators for both cases (purely continuous random variables and mixed random variables) that achieve the parametric MSE convergence rate regardless of the dimension as long as the densities are sufficiently smooth.

III-C Minimax Rate for MI estimation

We wrap up this section with a minimax lower bound on the MSE rate of convergence for the continuous MI estimation problem.

Theorem 4** (Bound on the Minimax Rate for Continuous $\mathbf{X},\mathbf{Y}$).**

Assume that $g$ is at least twice differentiable and that given $\epsilon>0,$ $\left|g^{\prime\prime}(\epsilon)\right|>0.$ Define the set of functions $\Sigma\left(s,H,\epsilon_{0},\epsilon_{\infty}\right)$ to be the set of Hölder continuous functions $\Sigma(s,H)$ that are bounded between $\epsilon_{0}$ and $\epsilon_{\infty}$. Then with $\gamma=\min\left\{8s/\left(4s+d_{X}+d_{Y}\right),1\right\}$, there exists a strictly positive constant $c$ such that

[TABLE]

The proof uses Le Cam’s method [79] and is given in Appendix F. Theorem 4 indicates that the minimax rate is the parametric rate $N^{-1}$ as long as $s\geq\left(d_{X}+d_{Y}\right)/4$. This is consistent with minimax rates for divergence [7, 8] and entropy [74] functional estimation, thus expanding the previous theory on minimax estimation of information theoretic functionals.

In Section V and Appendix B, we derive MI estimators that achieve the minimax rate when $s\geq d_{X}+d_{Y}$ and $s>\left(d_{X}+d_{Y}\right)/2$, respectively. While estimators have been derived for the divergence estimation problem that achieve the minimax rate for less smooth densities, they require numerical integration and are thus computationally slow [7, 8]. Deriving estimators of these functionals (e.g. MI and divergence) that are known to achieve the minimax rate in this less smooth regime and that are computationally reasonable thus remains an open problem.

IV Mixed Random Variables

In this section, we extend the results of Section III to general MI estimation when $\mathbf{X}$ and $\mathbf{Y}$ may have a mixture of discrete and continuous components. We focus on the most complex case: $\mathbf{X}$ and $\mathbf{Y}$ both have discrete and continuous components. The MI between $\mathbf{X}$ and $\mathbf{Y}$ is written in (4).

IV-A KDE Plug-in Estimator

We first define the KDE plug-in estimator of (4). Let $\mathcal{S}_{Y_{C}}$ and $\mathcal{S}_{X_{C}}$ be the respective supports of the corresponding densities of $\mathbf{Y}_{C}$ and $\mathbf{X}_{C}$ and let $\mathcal{S}_{Y_{D}}$ and $\mathcal{S}_{X_{D}}$ be the respective supports of the corresponding probability mass functions of $\mathbf{Y}_{D}$ and $\mathbf{X}_{D}$. Suppose we have $N$ i.i.d. samples of $(\mathbf{X},\mathbf{Y})$ drawn from $f_{XY}$ where the $i$th samples are denoted as $\left(\mathbf{X}_{i},\mathbf{Y}_{i}\right)=\left(\mathbf{X}_{i,C},\mathbf{X}_{i,D},\mathbf{Y}_{i,C},\mathbf{Y}_{i,D}\right)$. Define the following random variables:

[TABLE]

where $x\in\mathcal{S}_{X_{D}}$, $y\in\mathcal{S}_{Y_{D}}$, and $1_{\{\cdot\}}$ is the indicator function. These will be used to estimate the pmfs of the discrete components of $\mathbf{X}$ and $\mathbf{Y}$.

For the continuous components, we will condition on the discrete components and construct KDEs for the conditional probability density functions. Let $\mathcal{S}_{X_{C}}$ and $\mathcal{S}_{Y_{C}}$ be the respective supports of the marginal densities $f_{X_{C}}$ and $f_{Y_{C}}$ with corresponding dimensions of $d_{X}$ and $d_{Y}$. As before, let $K_{X}(\cdot)$ and $K_{Y}(\cdot)$ be kernel functions with $\int K_{X}(x)dx=\int K_{Y}(y)dy=1$, $||K_{X}||_{\infty},\,||K_{Y}||_{\infty}<\infty$ where $||K||_{\infty}=\sup_{x}|K(x)|$. Consider the following sets:

[TABLE]

The set $\mathcal{X}_{x}$ is the set of the continuous $\mathbf{X}$ data points where the corresponding discrete component is equal to $x$. The set $\mathcal{Y}_{y}$ is defined similarly. The KDEs for $f_{X_{C}|X_{D}}$, $f_{Y_{C}|Y_{D}}$, and $f_{X_{C}Y_{C}|X_{D}Y_{D}}$at $x\in\mathcal{S}_{X_{D}}$ and $y\in\mathcal{S}_{Y_{D}}$ are, respectively,

[TABLE]

where $\mathbf{Z}_{C}=\left(\mathbf{X}_{C},\mathbf{Y}_{C}\right)$ and $h_{Z_{C}|z}=\left(h_{X_{C}|x},h_{Y_{C}|y}\right)$. Note that we allow the bandwidths to depend on the discrete components of $\mathbf{X}$ and $\mathbf{Y}$. The reason for this is that the bandwidth is generally chosen as a function of the number of data points, which will differ for these conditional distributions as the discrete components of $\mathbf{X}$ and $\mathbf{Y}$ differ.

The MI $I(\mathbf{X};\mathbf{Y})$ can then be estimated by plugging in the conditional KDEs. Note that the MI in eq. (4) is written as a weighted sum of integral functionals. We therefore first define an intermediate estimator of the integral functionals:

[TABLE]

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

[TABLE]

The quality of the conditional density estimates in terms of bias and variance depends on the choice of bandwidths $h_{X_{C}|x}$ and $h_{Y_{C}|y}$. That is, for the KDE $\tilde{\mathbf{f}}_{X_{C}|x,h_{X_{C}|x}}$ to converge in MSE, it is necessary that $h_{X_{C}|x}\rightarrow 0$ and $\mathbf{N}_{x}h_{X_{C}|x}^{d_{X}}\rightarrow 0$ as $\mathbf{N}_{x}\rightarrow\infty$ (a similar result holds for $h_{Y_{C}|y}$) [80]. Furthermore, we will see when we derive the bias and variance of $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ that these conditions are also necessary for $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ to converge in MSE. Thus, when deriving the MSE convergence rate of $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$, we will assume that $h_{X_{C}|x}$ is a function of $\mathbf{N}_{x}$ and $h_{Y_{C}|y}$ is a function of $\mathbf{N}_{y}$.

IV-B MSE Convergence Rates of the Mixed Plug-in Estimator

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

Lemma 5**.**

Let $\mathbf{N}_{y}$, $\mathbf{N}_{x}$, and $\mathbf{N}_{xy}$ be defined as in (13). Assume that their corresponding probability mass functions are bounded away from zero. If $\alpha\in\mathbb{R}\backslash\{0,1\}$ and $\lambda+\beta+\gamma\in\mathbb{R}\backslash\{0,1\}$, then

[TABLE]

The proof is in Appendix G-A and uses the generalized binomial theorem, Taylor series expansions, and known results about the central moments of binomial random variables [81]. Lemma 5 provides key results on moments of products of the binomial random variables $\mathbf{N}_{xy}$, $\mathbf{N}_{x}$, and $\mathbf{N}_{y}$. These results can be used to derive the bias and variance of a plug-in estimator of MI with mixed components in (4) as long as the bias and variance of the corresponding plug-in estimator for the continuous weighted case in (6) is known. This is demonstrated in the following theorems for the KDE plug-in estimator $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$.

Theorem 6** (Bias Expansion for Mixed $\mathbf{X},\mathbf{Y}$).**

Assume that assumptions $\mathcal{A}.0$-$\mathcal{A}.5$ hold with respect to the functional $g$, the kernels $K_{X}$ and $K_{Y}$, and the densities $f_{X_{C}|X_{D}}$, $f_{Y_{C}|Y_{D}}$ and $f_{X_{C}Y_{C}|X_{D}Y_{D}}$. Assume that $\left|\mathcal{S}_{X_{D}}\right|,\left|\mathcal{S}_{Y_{D}}\right|<\infty$. Assume that $\mathbf{h}_{X_{C}|x}=l_{X}\mathbf{N}_{x}^{-\beta}$ and $\mathbf{h}_{Y_{C}|y}=l_{Y}\mathbf{N}_{y}^{-\alpha}$ with $0<\beta<\frac{1}{d_{X}}$, $0<\alpha<\frac{1}{d_{Y}}$, and $l_{X},l_{Y}>0$. Then the bias of $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ is

[TABLE]

The constants depend on the underlying densities, the chosen kernels, the functional $g$, and the probability mass functions and are independent of $l_{X}$, $l_{Y}$, and $N$. Furthermore, these rates are asymptotically tight.

Theorem 7** (Variance Bound for Mixed $\mathbf{X},\mathbf{Y}$).**

Assume that $\mathbf{h}_{X_{C}|x}=l_{X}\mathbf{N}_{x}^{-\beta}$ and $\mathbf{h}_{Y_{C}|y}=l_{Y}\mathbf{N}_{y}^{-\alpha}$ with $0<\beta<\frac{1}{d_{X}}$, $0<\alpha<\frac{1}{d_{Y}}$, $\beta d_{X}+\alpha d_{Y}\leq 1$, and $l_{X},l_{Y}>0$. Assume that $\left|\mathcal{S}_{X_{D}}\right|,\left|\mathcal{S}_{Y_{D}}\right|<\infty$. If the shaping function $g$ is Lipschitz continuous in both of its arguments, then the variance of $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ is $O(1/N)$.

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

IV-C Proof Sketches of Theorems 6 and 7

For Theorem 6, the proof splits the bias term into two terms by adding and subtracting $g\left(\mathcal{T}(\mathbf{X},\mathbf{Y})\frac{\mathbf{N}_{x}\mathbf{N}_{y}}{N\mathbf{N}_{xy}}\right)$ for each pair $(x,y)$ where $\mathcal{T}(\mathbf{X},\mathbf{Y})$ is independent of the data samples and is defined in Eq. (80). It can be shown that the newly added term has bias $O(1/N)$. The other term is handled by conditioning on the discrete components of the data samples to obtain the conditional bias terms $\mathbb{B}\left[\left.\tilde{\mathbf{G}}_{h_{X_{C}|x},h_{Y_{C}|y}}\right|\mathbf{X}_{1,D},\dots,\mathbf{X}_{N,D},\mathbf{Y}_{1,D},\dots,\mathbf{Y}_{N,D}\right]$ for each pair $(x,y)$. Theorem 2 can then be applied to each of these terms to obtain expressions of the random variables $\mathbf{N}_{x}$, $\mathbf{N}_{y}$, and $\mathbf{N}_{xy}$ with terms of the form given in Lemma 5. Lemma 5 can be applied to these terms to obtain the final result, where care is taken to ensure that all relevant terms have been handled properly. The full proof is given in Appendix G-B.

To prove Theorem 7, we use the law of total variance to split the variance into two terms: the expected value of the variance conditioned on the discrete components of the data samples and the variance of the conditional expectation. Theorem 3 is applied to the conditional variance term. For the conditional expectation term, we use results obtained in the proof of Theorem 6 combined with the Efron-Stein inequality [82] to obtain expressions of the random variables $\mathbf{N}_{x}$, $\mathbf{N}_{y}$, and $\mathbf{N}_{xy}$. Lemma 5 can be applied again to these terms to obtain the final result. The full proof is given in Appendix G-C.

V Ensemble Estimation of Generalized MI

If no bias correction is performed, then Theorems 2 and 6 show that the optimal bias rate of the KDE plug-in estimators $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ and $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ is $O\left(1/N^{1/(d_{X}+d_{Y}+1)}\right)$, which converges very slowly to zero when either $d_{X}$ or $d_{Y}$ are not small. Thus the standard KDE plug-in estimators will perform poorly in higher-dimensional settings. We use the theory of optimally weighted ensemble estimation developed in [14] to improve this rate. For brevity, we focus on the case where $\mathbf{X}$ and $\mathbf{Y}$ both contain a mixture of discrete and continuous components. The purely continuous case is described in Appendix B-A.

An ensemble of estimators is first formed by choosing different bandwidth values for the plug-in estimators as follows. Let $\mathcal{L}$ be a set of real positive numbers with $|\mathcal{L}|=L<\infty$. This set will parameterize the bandwidths $\mathbf{h}_{X_{C}|x}$ and $\mathbf{h}_{Y_{C}|y}$ for $\tilde{\mathbf{f}}_{X_{C}|x,h_{X_{C}|x}}$ and $\tilde{\mathbf{f}}_{Y_{C}|y,h_{Y_{C}|y}}$, respectively, resulting in $L$ estimators in the ensemble. In other words, we set $\mathbf{h}_{X_{C}|x}(l)=l\mathbf{N}_{x}^{-\beta}$ and $\mathbf{h}_{Y_{C}|y}(l)=l\mathbf{N}_{y}^{-\alpha}$ with $l\in\mathcal{L}$. While different parameter sets for $\mathbf{h}_{X_{C}|x}$ and $\mathbf{h}_{Y_{C}|y}$ can be chosen, we only use one set here for simplicity of exposition. To achieve the parametric rate, we need to ensure that the final terms in (26) are $O(1/\sqrt{N})$. Thus we require the following conditions to be met:

[TABLE]

For all of these conditions to hold, it is necessary that $s\geq d_{X}+d_{Y}$. Thus for each estimator in the ensemble we choose $\mathbf{h}_{X_{C}|x}(l)=l\mathbf{N}_{x}^{-1/\left(2\left(d_{X}+d_{Y}\right)\right)}$ and $\mathbf{h}_{Y_{C}|y}(l)=l\mathbf{N}_{y}^{-1/\left(2\left(d_{X}+d_{Y}\right)\right)}$ 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

[TABLE]

This is the weighted ensemble estimator. From Theorem 6, the bias of $\tilde{\mathbf{G}}_{w}$ is

[TABLE]

where we use $\theta$ notation to omit the constants.

We use the general theory of optimally weighted ensemble estimation in [73, 14] to improve the MSE convergence rate of the plug-in estimator by choosing the appropriate weights to cancel the lower order terms in (28):

Theorem 8** (Ensemble MSE).**

Let $\mathcal{L}$ be a set of real positive numbers with $|\mathcal{L}|=L<\infty$ and let $J=\{1,2,\dots,d_{X}+d_{Y}\}$. Assume the same conditions in Theorems 6 and 7 hold with $\mathbf{h}_{X_{C}|x}(l)=l\mathbf{N}_{x}^{-1/\left(2\left(d_{X}+d_{Y}\right)\right)}$ and $\mathbf{h}_{Y_{C}|y}(l)=l\mathbf{N}_{y}^{-1/\left(2\left(d_{X}+d_{Y}\right)\right)}$. Assume that $s\geq d_{X}+d_{Y}$ and define $\tilde{\mathbf{G}}_{w}$ as in (27). Then the MSE of $\tilde{\mathbf{G}}_{w_{0}}$ attains the parametric rate of convergence of $O\left(1/N\right)$ where $w_{0}$ is the solution to the following offline convex optimization problem:**

[TABLE]

To summarize, if the weights are chosen using eq. (29), then the weighted ensemble estimator $\tilde{\mathbf{G}}_{w_{0}}$ achieves the parametric MSE rate. In practice, the optimization problem in (29) typically results in a very large increase in variance for finite samples. Thus we use a relaxed version of (29):

[TABLE]

The parameter $\eta$ is chosen to achieve a trade-off between bias and variance. As shown in [16, 14], the ensemble estimator $\tilde{\mathbf{G}}_{w_{0}}$ using the resulting weight vector from the optimization problem in (30) still achieves the parametric MSE convergence rate under the same assumptions as described previously. We denote this estimator as $\tilde{\mathbf{G}}_{GENIE}$. Algorithm 1 summarizes the estimator $\tilde{\mathbf{G}}_{GENIE}$.

A similar approach can be used to derive an ensemble estimator for the case when $\mathbf{X}$ and $\mathbf{Y}$ are purely continuous. Furthermore, we can define ensemble estimators for both the continuous and the mixed cases that achieve the parametric MSE rate if $s>\left(d_{X}+d_{Y}\right)/2$, although the optimization problem is more complicated. See Appendix B for details.

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

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

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

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

V-A Parameter Selection

Asymptotically, 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.

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. 2.

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

Select $L=|\mathcal{L}|>|J|=I$ to obtain a feasible solution for the optimization problems in (29) and (30). We find that choosing a value of $30\leq L\leq 60$, 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. [14] gives more details on ensemble estimator parameter selection for continuous divergence estimation. These details also apply to the continuous parts of the mixed cases for MI estimation in this paper. In particular, the minimum and maximum bandwidth parameters can be efficiently selected based on the $k$ nearest neighbor distances of all data points.

Since the optimal weight $w_{0}$ can be calculated offline, the computational complexity of the estimators is dominated by the construction of the KDEs which has a complexity of $O\left(N^{2}\right)$ using the standard implementation. For very large datasets, more efficient KDE implementations (e.g. [85]) can be used to reduce the computational burden.

V-B Central Limit Theorem

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

Theorem 9** (CLT for Continuous $\mathbf{X},\mathbf{Y}$).**

Let $\tilde{\mathbf{G}}_{w}^{cont}$ be a weighted KDE ensemble estimator of $I_{\nu}(\mathbf{X};\mathbf{Y})$ when $\mathbf{X}$ and $\mathbf{Y}$ are continuous with bandwidths $h_{X}(l)$ and $h_{Y}(l)$ for each estimator in the ensemble. Assume that the shaping function $g$ is Lipschitz in both arguments with Lipschitz constant $C_{g}$ and that $h_{X}(l),\,h_{Y}(l)\rightarrow 0$, $N\rightarrow\infty$, and $Nh_{X}^{d_{X}}(l),\,Nh_{Y}^{d_{Y}}(l)\rightarrow\infty$ for each $l\in\mathcal{L}$. Then for fixed $\mathcal{L}$, and if $\mathbf{S}$ is a standard normal random variable,

[TABLE]

The proof is based on an application of Slutsky’s Theorem preceded by an application of the Efron-Stein inequality (see Appendix H).

For the mixed component case, if $\mathcal{S}_{X}$ and $\mathcal{S}_{Y}$ are finite, then the corresponding ensemble estimators also obey a central limit theorem. The proof follows by an application of Slutsky’s Theorem combined with Theorem 9.

Corollary 10** (CLT for mixed $\mathbf{X},\mathbf{Y}$).**

Let $\tilde{\mathbf{G}}_{w}$ be a weighted KDE ensemble estimator of $I(\mathbf{X};\mathbf{Y})$ when $\mathbf{X}$ and $\mathbf{Y}$ contain both continuous and discrete components. Let the bandwidths for the conditional estimators be $\mathbf{h}_{X_{C}|x}(l)$ and $\mathbf{h}_{Y_{C}|y}(l)$ for each estimator in the ensemble. Assume that the shaping function $g$ is Lipschitz in both arguments and that $\mathbf{h}_{X_{C}|x},\,\mathbf{h}_{Y_{C}|y}\rightarrow 0$, $N\rightarrow\infty$, and $Nh_{X}^{d_{X}},\,Nh_{X|y}^{d_{X}}\rightarrow\infty$ for each $l\in\mathcal{L}$ and $\forall(x,y)\in\mathcal{S}_{X_{D}}\times\mathcal{S}_{Y_{D}}$ with $\left|\mathcal{S}_{X_{D}}\right|,\left|\mathcal{S}_{Y_{D}}\right|<\infty$. Then for fixed $\mathcal{L}$,

[TABLE]

VI Applications

VI-A Simulations

In this section, we validate our theory by estimating the Rényi-$\alpha$ MI integral (i.e. $g(x)=x^{\alpha}$ in (3); see [49]) where $\mathbf{X}$ is a mixture of truncated Gaussian random variables restricted to the unit cube and $\mathbf{Y}$ is a categorical random variable that indicates the corresponding truncated Gaussian random variable that $\mathbf{X}$ is drawn from in the mixture. In this setting, $\mathbf{Y}$ can be viewed as a classification variable and $\mathbf{X}$ contains the chosen features, which are all continuous in this case. Since $\mathbf{X}$ is purely continuous and $\mathbf{Y}$ is purely discrete, the MI integral reduces to the following:

[TABLE]

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

We consider two cases. In the first case, $\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)}$.

Figure 1 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 $\tilde{\mathbf{G}}_{GENIE}$ for various sample sizes and for $d=4,\,6,\,9$, respectively. The ensemble estimator GENIE outperforms the standard plug-in estimator, especially for larger sample sizes and larger dimensions. This demonstrates that while an individual kernel estimator performs poorly, an ensemble of estimators including the individual estimator performs well.

For the second case, $\mathbf{Y}$ has six possible outcomes (i.e. $|\mathcal{S}_{Y}|=6$) and respective probabilities $\Pr(\mathbf{Y}=0)=0.35$, $\Pr(\mathbf{Y}=1)=0.2$, $\Pr(\mathbf{Y}=2)=\Pr(\mathbf{Y}=3)=0.15$, $\Pr(\mathbf{Y}=4)=0.1$, and $\Pr(\mathbf{Y}=5)=0.05$. We chose $\alpha=0.5$ and $d=6$. The conditional covariance 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 results are again given in Figure 1. The parameters for the ensemble estimator and the KDE plug-in estimators are the same as in the other three plots in Figure 1. The ensemble estimator also outperforms the plug-in estimator in this setting.

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

VI-B Application to Single-Cell RNA-Sequencing Data

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

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

For these datasets, we estimated two MI measures: the Rényi MI and DREMI [37]. We define the Rényi MI to be equal to the Rényi divergence between the joint distribution of $\mathbf{X}$ and $\mathbf{Y}$ and the product of the marginal distributions. The DREMI score is a weighted MI developed specifically for analyzing single-cell data [37]. See Appendix A-A for further details. Note that no other estimator has been defined for $I_{DREMI}$ when the dimension of the continuous component or components are greater than 1.

VI-B1 Mouse bone marrow data

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

For this data, we estimated the two different MI measures between the cell type classification (discrete) and selected groups of genes (continuous). We estimated the MI for different combinations of genes selected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with the hematopoietic cell lineage [95, 96, 97]. Each of these collections contained 8-10 genes. Since the number of cell types is discrete and the gene expression levels are continuous, the estimation problem corresponds to estimating the MI between $\mathbf{X}$ and $\mathbf{Y}$ for the case where $\mathbf{Y}$ is discrete and $\mathbf{X}$ is continuous. In this problem, $|\mathcal{S}_{Y}|=19$ and $d_{X}$ is the number of genes in the chosen collection.

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

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

VI-B2 Human embryoid body data

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

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

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

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

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

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

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

VII Conclusion

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

Appendix A Extended Genomics Details and Results

Here we provide further details on the genomics experiments.

A-A Imputation and DREMI

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

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

[TABLE]

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

A-B EB Data Extended Results

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

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

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

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

Appendix B MI Ensemble Estimation Extensions

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

B-A Continuous Random Varables

We can apply Theorem 3 in [14] to obtain a version of the GENIE MI estimator that achieves the parametric rate for the case when $\mathbf{X}$ and $\mathbf{Y}$ are purely continuous. For completeness, we repeat the theorem and its proof here. For a general estimation problem, let $N$ be the number of available samples and let $\mathcal{L}=\left\{l_{1},\dots,l_{L}\right\}$ be a set of index values. For an indexed ensemble of estimators $\left\{\hat{\mathbf{E}}_{l}\right\}_{l\in\mathcal{L}}$ of a parameter $E$, the weighted ensemble estimator with weights $w=\left\{w\left(l_{1}\right),\dots,w\left(l_{L}\right)\right\}$ satisfying $\sum_{l\in\mathcal{L}}w(l)=1$ is defined as

[TABLE]

Consider the following conditions on $\left\{\hat{\mathbf{E}}_{l}\right\}_{l\in\mathcal{L}}$:

• •

$\mathcal{C}.1$ The bias is expressible as

[TABLE]

where $c_{i}$ are constants depending on the underlying density and are independent of $N$ and $l$, $J=\left\{i_{1},\dots,i_{I}\right\}$ is a finite index set with $I<L$, and $\psi_{i}(l)$ are basis functions depending only on the parameter $l$ and not on the sample size $N$.

• •

$\mathcal{C}.2$ The variance is expressible as

[TABLE]

Theorem 11** (Theorem 3 in [14]).**

Assume 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 the weighted ensemble estimator attains the parametric rate of convergence:

[TABLE]

The weight vector $w_{0}$ is the solution to the following convex optimization problem:**

[TABLE]

Proof.

Due to condition $\mathcal{C}.1$, the bias of the weighted ensemble estimator is

[TABLE]

Denote the covariance matrix of the ensemble of estimators as $\Sigma_{L}$. By the Cauchy-Schwarz inequality and condition $\mathcal{C}.2$, the entries of $\Sigma_{L}$ are $O(1/N)$. The variance of the weighted estimator is then bounded above as

[TABLE]

The optimization problem in (31) zeroes out the lower-order bias terms and limits the $\ell_{2}$ norm of the weight vector $w$ to prevent the variance from exploding. This results in an MSE of $O(1/N)$ as long as the dimension $d$ is fixed and $L$ is fixed and independent of the sample size $N$. A solution is guaranteed to (31) as long as $L>I$ and the vectors $a_{i}=[\psi_{i}(l_{1}),\dots,\psi_{i}(l_{L})]$ are linearly independent.

∎

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

[TABLE]

We can use (34) to obtain a GENIE estimator for the purely continuous case. Theorem 2 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 2 and 3, 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 (34). The resulting estimator

[TABLE]

achieves the parametric MSE rate when $s\geq d_{X}+d_{Y}$. We denote this estimator as $\tilde{\mathbf{G}}_{GENIE}^{cont}$.

B-B Less Smooth Densities

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

B-B1 Continuous Random Variables

We first consider the case where $\mathbf{X}$ and $\mathbf{Y}$ are both purely continuous. Consider the following result on the bias of the plug-in estimator:

Theorem 12** (Continuous Bias Expanded).**

Assume that assumptions $\mathcal{A}.0$-$\mathcal{A}.5$ hold. Then for $\lambda\geq 2$ a positive integer, the bias of $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ is

[TABLE]

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

We now use these results to define a new ensemble estimator. Set $\delta>0$ and 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/(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\geq 1$ and $i+j+m\delta\leq\frac{d_{X}+d_{Y}+\delta}{2}.$ These all correspond to terms that converge to zero slower than $N^{-1/2}$ when left uncorrected. The optimal weight $w_{0}$ is again calculated using (34) and the resulting ensemble estimator 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$.

B-B2 Mixed Random Variables

We now consider the case where $\mathbf{X}$ and $\mathbf{Y}$ may have any mixture of continuous and discrete components. We have a similar result on the bias as in Theorem 12. Here we assume that $\mathbf{h}_{X_{C}|x}=l_{X}\mathbf{N}_{x}^{-\beta}$ and $\mathbf{h}_{Y_{C}|y}=l_{Y}\mathbf{N}_{y}^{-\alpha}$ with $0<\beta<\frac{1}{d_{X}}$, $0<\alpha<\frac{1}{d_{Y}}$, and $l_{X},l_{Y}>0$.

Theorem 13** (Mixed Bias Expanded).**

Assume that the same assumptions hold as in Theorem 6. Then for $\lambda\geq 2$ a positive integer, the bias of $\tilde{\mathbf{G}}_{h_{X_{C}|X_{D}},h_{Y_{C}|Y_{D}}}$ is

[TABLE]

The proof is given in Appendix G-B.

We now use these results to define a new ensemble estimator in the mixed case. The procedure is similar to the continuous case. Set $\delta>0$ and 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 $\mathbf{h}_{X_{C}|x}(l_{X})=l_{X}\mathbf{N}_{x}^{-1/(d_{X}+d_{Y}+\delta)}$ and $\mathbf{h}_{Y_{C}|y}(l_{Y})=l_{Y}\mathbf{N}_{y}^{-1/(d_{X}+d_{Y}+\delta)}$. 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$. The first set of terms in the optimization problem are $\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}$. The second set of terms are $\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\geq 1$ and $i+j+m\delta\leq\frac{d_{X}+d_{Y}+\delta}{2}.$ The optimal weight $w_{0}$ is again calculated using (34) and the resulting ensemble estimator 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$.

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

Appendix C The Boundary Condition

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

Theorem 14**.**

Let $K$ be a kernel function with rectangular or circular support, i.e., $K(u)=0$ for $\left\|u\right\|_{1}>1$ or $\left\|u\right\|_{2}>1$, respectively. Let $K\in\Sigma(s,H_{K})$ in the interior of its support. 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. Then for $\mathcal{S}=[0,1]^{d}$ and any positive integers $t$ and $m$, we have that

[TABLE]

where $v_{t}(h)$ admits the expansion

[TABLE]

for some constants $e_{i,q,t,l}$.

Before proving this, we will relate this result to assumption $\mathcal{A}.5$.

Corollary 15**.**

Assumption $\mathcal{A}.5$ is satisfied under the same conditions as in Theorem 14, i.e., $K\in\Sigma(s,H_{K})$ is a kernel function with either rectangular or circular support and $\mathcal{S}=[0,1]^{d}$.

This follows immediately from assumption $\mathcal{A}.2$ that $f$ is in the Hölder class $\Sigma(s,H)$, which implies that $D^{q}f(x)$ is $r-|q|$ times differentiable. Thus equation (37) implies $\mathcal{A}.5$.

C-A Proof of Theorem 14: Rectangular Support Kernels

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

Single Coordinate Boundary Point

Since $K\in\Sigma(s,H_{K})$, we can take a Taylor series expansion of $K^{l}(u)$ around zero to obtain

[TABLE]

where $p_{K,l}(u)$ is a polynomial function of $u$ with degree $r$.

Consider points $x$ that are boundary points by virtue of a single coordinate $x_{i}$ such that $x_{i}+u_{i}h\notin\mathcal{S}$. Without loss of generality, assume that $x_{i}+u_{i}h>1$. After performing the above Taylor series expansion of $K^{l}$, the inner integral in (37) can then be evaluated first with respect to all coordinates other than $i$. Since all of these coordinates lie within the support, the inner integral over these coordinates will amount to integration of the polynomial $p_{x}(u)$ over a symmetric $d-1$ dimensional rectangular region $\left|u_{j}\right|\leq 1$ for all $j\neq i$. This yields a function $\sum_{m=1}^{q+r}\tilde{p}_{m}(x)u_{i}^{m}+o(u_{i}^{r+1})$ where the coefficients $\tilde{p}_{m}(x)$ are each $r-|q|$ times differentiable wrt $x$.

With respect to the $u_{i}$ coordinate, the inner integral will have limits from $\frac{1-x_{i}}{h}$ to $1$ for some $1>x_{i}>1-h$. Consider the $\tilde{p}_{m}(x)u_{i}^{m}$ monomial term. The inner integral wrt this term yields (ignoring the $o(\cdot)$ term for now)

[TABLE]

Raising the right hand side of (39) to the power of $t$ results in an expression of the form

[TABLE]

where the coefficients $\check{p}_{j}(x)$ are $r-|q|$ times differentiable wrt $x$. Integrating (40) over all the coordinates in $x$ other than $x_{i}$ results in an expression of the form

[TABLE]

where again the coefficients $\bar{p}_{j}(x_{i})$ are $r-|q|$ times differentiable wrt $x_{i}$. Note that since the other cooordinates of $x$ other than $x_{i}$ are far away from the boundary, the coefficients $\bar{p}_{j}(x_{i})$ are independent of $h$. To evaluate the integral of (41), consider the $r-|q|$ term Taylor series expansion of $\bar{p}_{j}(x_{i})$ around $x_{i}=1$, where we use a smooth extension of the function and its derivatives to the boundary. This will yield terms of the form

[TABLE]

for $0\leq j\leq r-|q|$, and $0\leq k\leq(|q|+r)t$. By a similar analysis, the $o(\cdot)$ terms from the Taylor expansion of the kernel result in $o(h^{r})$. Combining terms results in the expansion $v_{t}(h)=\sum_{i=1}^{r-|q|}e_{i,q,t}h^{i}+o\left(h^{r-|q|}\right)$.

Multiple Coordinate Boundary Point

The case where multiple coordinates of the point $x$ are near the boundary is a straightforward extension of the single boundary point case so we only sketch the main ideas here. As an example, consider the case where 2 of the coordinates are near the boundary. Assume for notational ease that they are $x_{1}$ and $x_{2}$ and that $x_{1}+u_{1}h>1$ and $x_{2}+u_{2}h>1$. The inner integral in (37) can again be evaluated first wrt all coordinates other than 1 and 2. This yields a function $\sum_{m,j=1}^{q+r}\tilde{p}_{m,j}(x)u_{1}^{m}u_{2}^{j}$ where the coefficients $\tilde{p}_{m,j}(x)$ are each $r-q$ times differentiable wrt $x$. Integrating this wrt $x_{1}$ and $x_{2}$ and then raising the result to the power of $t$ yields a double sum similar to (40). Integrating this over all the coordinates in $x$ other than $x_{1}$ and $x_{2}$ gives a double sum similar to (41). Then a Taylor series expansion of the coefficients and integration over $x_{1}$ and $x_{2}$ yields the result.

C-B Proof of Theorem 14: Circular Support Kernels

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

Single Coordinate Boundary Point

Consider points $x$ that are boundary points due to a single coordinate $x_{i}$ s.t. $x_{i}+u_{i}h\notin\mathcal{S}$. Without loss of generality, assume that $x_{i}+u_{i}h>1$. We focus first on the inner integral in (37). We will use the following lemma:

Lemma 16**.**

Let $D_{d}(\rho)$ be a $d$-sphere with radius $\rho$ and let $\sum_{i=1}^{d}n_{i}=q$. Then

[TABLE]

where $C$ is a constant that depends on the $n_{i}$s and $d$.

Proof.

We convert to $d$-dimensional spherical coordinates to handle the integration. Let $r$ be the distance of a point $u$ from the origin. We nave $d-1$ angular coordinates $\phi_{i}$ where $\phi_{d-1}$ ranges from [math] to $2\pi$ and all other $\phi_{i}$ range from [math] to $\pi$. The conversion from the spherical coordinates to Cartesian coordinates is then

[TABLE]

The spherical volume element is then

[TABLE]

Combining these results gives

[TABLE]

∎

Since $K\in\Sigma(s,H_{K})$, we can take a Taylor series expansion of $K^{l}(u)$ around zero to obtain

[TABLE]

where $p_{K,l}(u)$ is a polynomial function of $u$ with degree $r$.

The region of integration for the inner integral in (37) corresponds to a hyperspherical cap with radius $1$ and height of $\frac{1-x_{i}}{h}$. The inner integral can be calculated using an approach similar to that used in [123] to calculate the volume of a hyperspherical cap. It is obtained by integrating the polynomial $p_{x}(u)$ over a $d-1$-sphere with radius $\sin\theta$ and height element $d\cos\theta$. This is done using Lemma 16. We then integrate over $\theta$ which has a range of [math] to $\phi=\cos^{-1}\left(\frac{1-x_{i}}{h}\right).$ Thus we have

[TABLE]

where $\tilde{p}_{m}(x)$ is the polynomial coefficient corresponding to $u_{i}^{m}$ after the Taylor series expansion of $K^{l}$ and integrating over the $d-1$-sphere.

We will focus on the first term in Eq. 45 as the $o(\cdot)$ term will follow similarly as a polynomial function of $u$. From standard integral tables, we get that for $n\geq 2$ and $m\geq 0$

[TABLE]

If $n=1$, then we get

[TABLE]

Since $\phi=\cos^{-1}\left(\frac{1-x_{i}}{h}\right)$, we have

[TABLE]

Therefore, if $n$ is odd, we obtain

[TABLE]

where the constants depend on $m$ and $n$.

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

[TABLE]

If $m=2$, then

[TABLE]

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

[TABLE]

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

[TABLE]

If $d$ is odd, then combining (47) and (49) with (45) gives

[TABLE]

where the coefficients $p_{m,\ell,l}(x)$ are $r-|q|$ times differentiable wrt $x$. Similarly, if $d$ is even, then

[TABLE]

where again the coefficients $p_{m,\ell,l}^{{}^{\prime}}(x)$ and $p^{{}^{\prime}}(x)$ are $r-q$ times differentiable wrt $x$. Raising (50) and (51) to the power of $t$ gives respective expressions of the form

[TABLE]

[TABLE]

where the coefficients $\check{p}_{m,\ell,l}(x)$ and $\check{p}_{m,\ell,l,n}(x)$ are all $r-q$ times differentiable wrt $x$. Integrating (52) and (53) over all the coordinates in $x$ except for $x_{i}$ affects only the $\check{p}_{m,\ell,l}(x)$ and $\check{p}_{m,\ell,l,n}(x)$ coefficients, resulting in respective expressions of the form

[TABLE]

[TABLE]

The coefficients $\bar{p}_{m,\ell,l}(x_{i})$ and $\bar{p}_{m,\ell,l,n}(x_{i})$ are $r-q$ times differentiable wrt $x_{i}$. Since the other coordinates of $x$ other than $x_{i}$ are far away from the boundary, the coefficients are independent of $h$. For the integral wrt $x_{i}$ of (54), taking a Taylor series expansion of $\bar{p}_{m,\ell}(x_{i})$ around $x_{i}=1$ (again using a smooth extension of the function and its derivatives to the boundary) yields terms of the form

[TABLE]

where $0\leq j\leq r-q$, $0\leq\ell\leq(d+|q|)t$, $0\leq m\leq|q|t$, and $B(x,y)$ is the beta function. Note that the first step uses the substitution of $y_{i}=\left(\frac{1-x_{i}}{h}\right)^{2}$.

If $d$ is even (i.e. (55)), a simple closed-form expression is not easy to obtain due to the $\cos^{-1}\left(\frac{1-x_{i}}{h}\right)$ terms. However, by similarly applying a Taylor series expansion to $\bar{p}_{m,\ell,l,n}(x_{i})$ and substituting $y_{i}=\frac{1-x_{i}}{h}$ gives terms of the form of

[TABLE]

for $0\leq j\leq r-|q|$, $0\leq\ell\leq(d+|q|)t$, $0\leq m\leq(|q|+r)t$, and $0\leq n\leq t$. By a similar analysis, the $o(\cdot)$ term in Eq. 45 results in $o(h^{r})$. Combining terms results in the expansion $v_{t}(h)=\sum_{i=1}^{r-|q|}e_{i,q,t,l}h^{i}+o(h^{r-|q|})$.

Multiple Coordinate Boundary Point

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

C-C Proof of Proposition 1

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

Lemma 1**.**

[10]** Let the density support set $\mathcal{S}$ be the unit cube with $d\geq 2$ and let the derivatives of the density $f$ vanish near the boundary of the support set, which is denoted as $\partial\mathcal{S}$. Assume that $f$ belongs to the Hölder set of order $s$, i.e. $f\in\Sigma(H,s)$. Assume that $q$ with $0\leq|q|\leq r$ and $x\in\mathcal{B}:=\{x\in\mathcal{S}|\text{dist}(x,\partial\mathcal{S})\leq h\}$. Then if $h\leq 1/2$,

[TABLE]

Proof.

We include a proof here for completeness. We will use an induction argument on $|q|$ as $|q|$ decreases from $r$ to 0. Let $y\in\partial\mathcal{S}$ such that $||x-y||\leq h$. For $|q|=r$, we have

[TABLE]

where the first step comes from the fact that $D^{q}f(y)=0$, the second step uses the Hölder condition, and the last step follows from the choice of $y$. Suppose we have the desired bound for derivatives of order $|q|+1$. Let $x\in\mathcal{B}$ and let $u=(0,\dots,0,\pm 1,0,\dots,0)\in\mathbb{R}^{d}$, where $u_{j}=\pm 1$ for some $j\in[d]$. Then we can write $x=y+hu$ for some $y\in\mathcal{B}$. Furthermore, the point $y+tu\in\mathcal{B}$ for $t\in[0,h]$. To see how this is true, note that if $x\in\mathcal{B}$ then there exists at least one coordinate $x_{i}$ such that either $0\leq x_{i}\leq h$ or $1-h\leq x_{i}\leq 1$. Pick any other coordinate $x_{j}$. Then choose the sign of $u_{j}$ such that $y=x-hu\in\mathcal{S}$. This is possible since $h\leq 1/2$. Since $y$ is constructed by moving from $x\in\mathcal{B}$ in parallel with the boundary, then $y\in\mathcal{B}$ as well as $y+tu\in\mathcal{B}$.

Based on this construction of $y$ and $u$, we have

[TABLE]

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

∎

From Lemma 1,

[TABLE]

Consider the case when $K$ has rectangular support from (without loss of generality) $-1$ to 1 in each dimension. We consider the setting when $x$ is close to the 1 boundary in a single coordinate $x_{i}$. The inner integral in (57) then reduces to

[TABLE]

where $C$ is a constant that comes from integrating the other terms of $u^{q}$ and $m\leq|q|$ is the exponent of $u_{i}$ in $u^{q}$. Based on equation (43), taking the integral of this result with respect to $x$ will then yield terms of the form of $h^{j+1}$ for $j\geq 0$. The extension to the case with multiple boundary coordinates is similar to that in the proof of Theorem 14.

Combining this with the results in equations (57) and (59) gives that

[TABLE]

Therefore, assumption $\mathcal{A}.5$ is satisfied with the coefficients in the expansion being equal to zero.

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

C-D Assumption $\mathcal{A}.5$ and the Truncated Gaussian Distribution

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

[TABLE]

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

[TABLE]

Therefore, the integral with respect to $x$ in (62) has terms of the form of $\frac{(1-x)^{j+k}}{h^{k}}$ which integrates to $\frac{h^{j+1}}{j+k+1}$ (see (43)). Combining these results gives the desired expansion. As an example, suppose $t=1$ and $q=1$. Then (62) becomes

[TABLE]

Distributing the expression for $\frac{d}{dx}f(x)$ through and evaluating the integrals separately gives

[TABLE]

[TABLE]

Taking the difference between these expressions gives the final expansion:

[TABLE]

Appendix D Proof of Theorem 2 (Continuous Bias)

Here we prove the results shown in Theorem 12, which includes the results for Theorem 2. In this section and throughout all other proofs, let $\mathbb{E}_{\mathbf{Z}}$ denote the conditional expectation given $\mathbf{Z}$.

The bias of $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ is

[TABLE]

where $\mathbf{X}$ and $\mathbf{Y}$ are drawn jointly from $f_{XY}$. We will derive an expression for this in terms of the bandwidths by applying Taylor series expansions to both the functional $g$ and the densities.

The Taylor series expansion of $g\left(\frac{\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y})\nu_{1}\nu_{2}}{\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{X},\mathbf{Y})\nu_{3}}\right)$ around $f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})\nu_{1}\nu_{2}$ and $f_{XY}(\mathbf{X},\mathbf{Y})\nu_{3}$ gives an expansion with terms of the form of

Define $g\left(t_{1}/t_{2}\right)=g\left(t_{1},t_{2}\right).$ Define the following terms:

[TABLE]

Then the Taylor series expansion of $g\left(\frac{\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y})\nu_{1}\nu_{2}}{\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{X},\mathbf{Y})\nu_{3}}\right)$ around $f_{X}(\mathbf{X})f_{Y}(\mathbf{Y})\nu_{1}\nu_{2}$ and $f_{XY}(\mathbf{X},\mathbf{Y})\nu_{3}$ gives

[TABLE]

To simplify this, we focus on the terms in (67). Note that if $\nu_{i}=1$, then the terms in (67) are unaffected. For other values, $\nu_{i}^{j}$ decreases to zero as $j\rightarrow\infty$ since $0<\nu_{i}<1$.

By the binomial theorem, we obtain

[TABLE]

To derive the bias, we will take the conditional expectation given $\mathbf{Z}$ of these terms.

Since we are not doing explicit 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^{{}^{\prime}}\notin\mathcal{S}_{X}\times\mathcal{S}_{Y}$, $K_{X}\left(\frac{X-X^{{}^{\prime}}}{h_{X}}\right)K_{Y}\left(\frac{Y-Y^{{}^{\prime}}}{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.

Interior Points

We first consider the case where $\mathbf{Z}=(\mathbf{X},\mathbf{Y})$ is drawn from $f_{XY}$ in the interior of $\mathcal{S}_{X}\times\mathcal{S}_{Y}$. It can be shown (see [14]) by Taylor series expansions of the probability densities that

[TABLE]

The constants in the above expressions are independent of the bandwidths $h_{X}$ and $h_{Y}$ and only depend on the densities and their derivatives.

We will also need an expression for $\mathbb{E}_{\mathbf{Z}}\left[\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y})\right]$:

[TABLE]

where we have used the fact that $\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$ are independent when $i\neq j$.

We will also need expressions for the conditional expectation of powers of the KDEs to simplify (68). Consider first $\left(\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{Z})\right)^{2}$. Note that

[TABLE]

The above double sum can be split into two cases: $i=j$ ($N$ terms) and $i\neq j$ ($N^{2}-N$ terms). When $i=j$, we have by Taylor series expansions of the joint density $f_{XY}$:

[TABLE]

When $i\neq j,$recall that $\mathbf{Z}_{i}$ is independent of $\mathbf{Z}_{j}$. Therefore, we obtain for these terms

[TABLE]

Combining these results gives

[TABLE]

For the cross term ($k=1$) in (68) when $q=2$, we obtain

[TABLE]

Combining these results gives

[TABLE]

By following similar procedures, it can be shown that for $q\geq 2$

[TABLE]

In particular, the terms related to $f_{XY}(\mathbf{Z})^{q}$ all combine to be $O(1/N)$. In the $q=2$ example above, we end up with

[TABLE]

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

[TABLE]

A similar pattern holds for higher values of $q$.

Then by following a similar process, we obtain for $q\geq 2$

[TABLE]

Points Near the Boundary

For a point near the boundary of the support, we extend the expectation beyond the support of the density. As an example if $\mathbf{X}$ is near the boundary of $\mathcal{S}_{X}$, then we get

[TABLE]

Note that we are technically evaluating the density $f_{X}$ at points outside of the support in $T_{1,X}(\mathbf{X})$. However, to obtain an expression for this integral, we take a Taylor series expansion of $f_{X}$ at the point $\mathbf{X}$ which is inside the support. Thus the exact manner in which we define the extension of $f_{X}$ does not matter as long as the Taylor series remains the same and as long as the extension is smooth. Thus the expected value of $T_{1,X}(\mathbf{X})$ gives an expression similar to that of the interior point case in (69).

For the $T_{2,X}(\mathbf{X})$ term, we can use multi-index notation on the expansion of $f_{X}$ to show that

[TABLE]

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

[TABLE]

Similar expressions can be obtained for $\tilde{\mathbf{f}}_{Y,h_{Y}}$, $\tilde{\mathbf{f}}_{Z,h_{Z}}$, and the product $\tilde{\mathbf{f}}_{Y,h_{Y}}\tilde{\mathbf{f}}_{X,h_{X}}$.

The above results considers $\tilde{\mathbf{e}}_{Z,h_{Z}}^{q}(\mathbf{Z})$ and $\tilde{\mathbf{e}}_{XY,h_{X},h_{Y}}^{q}(\mathbf{Z})$ for $q=1$. We now consider when $q\geq 2$. We follow a similar procedure where we extend the density beyond the support, but only evaluate the densities and their derivatives at points within the support. Thus by the binomial theorem, we can write

[TABLE]

As before, $T_{1,q,Z}(\mathbf{Z})$ corresponds to the case where we have extended the density beyond the support and results in terms of the form in (71). $T_{2,q,Z}(\mathbf{Z})$ corresponds to the case where we integrate outside of the boundary. The additional powers applied to the KDE simply result in terms with the kernel raised to a power or $\mathbb{E}_{\mathbf{Z}}\left[\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{Z})\right]$ raised to a power. By applying assumption $\mathcal{A}.5$, we obtain $\mathbb{E}\left[T_{2,q,X}(\mathbf{X})\right]=\sum_{i=1}^{r}\sum_{j=1}^{r}e_{q,i}h_{X}^{i}h_{Y}^{j}+o\left(h_{X}^{r}+h_{Y}^{r}\right)$. Similar results are obtained for $\tilde{\mathbf{e}}_{XY,h_{X},h_{Y}}^{q}(\mathbf{Z})$.

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

Appendix E Proof of Theorem 3 (Continuous Variance)

Here we prove Theorem 3. The proof uses the Efron-Stein inequality [124, 82]:

Lemma 17** (Efron-Stein Inequality).**

Let $\mathbf{X}_{1},\dots,\mathbf{X}_{n},\mathbf{X}_{1}^{{}^{\prime}},\dots,\mathbf{X}_{n}^{{}^{\prime}}$ 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

[TABLE]

In this case we consider the samples $\left\{\mathbf{Z}_{1},\dots,\mathbf{Z}_{N}\right\}$ and $\left\{\mathbf{Z}_{1}^{{}^{\prime}},\mathbf{Z}_{2}\dots,\mathbf{Z}_{N}\right\}$ and the respective estimators $\tilde{\mathbf{G}}_{h_{X},h_{Y}}$ and $\tilde{\mathbf{G}}_{h_{X},h_{Y}}^{{}^{\prime}}$. By the triangle inequality,

[TABLE]

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

[TABLE]

[TABLE]

For the $\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{Z}_{1})$ term, we first apply Jensen’s inequality:

[TABLE]

By making the substitutions $\mathbf{u}_{i}=\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}$, $\mathbf{v}_{i}=\frac{\mathbf{Y}_{1}-\mathbf{Y}_{i}}{h_{Y}}$, $\mathbf{u}_{i}^{{}^{\prime}}=\frac{\mathbf{X}_{1}-\mathbf{X}_{i}}{h_{X}}$, and $\mathbf{v}_{i}^{{}^{\prime}}=\frac{\mathbf{Y}_{1}^{{}^{\prime}}-\mathbf{Y}_{i}}{h_{Y}}$ in the expectation, we obtain

[TABLE]

[TABLE]

This gives

[TABLE]

where we have used the fact that $\nu_{3}\leq 1$.

For the product of the marginal KDEs, we have that

[TABLE]

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

[TABLE]

For the second term in (73), it can be shown that (see [14])

[TABLE]

By a similar approach,

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

Applying Jensen’s inequality in conjunction with these results gives

[TABLE]

Applying the Efron-Stein inequality finishes the proof.

Appendix F Proof of Minimax Rates (Theorem 4)

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

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

[TABLE]

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

[TABLE]

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

Lemma 18**.**

Let $I$ be a functional defined on some subset of a parameter space $\Theta\times\Theta$ which contains $\left(f_{XY},f_{X}f_{Y}\right)$ and $\left(f_{XY,\lambda},f_{X}f_{Y}\right)$ for all $\lambda$ in some index set $\Lambda$. Denote the distributions of $f_{XY},f_{X}f_{Y},$ and $f_{XY,\lambda}$ as $F_{XY}$, $F_{X}F_{Y}$, and $F_{XY,\lambda}$, respectively. Define $\bar{F}_{XY}^{N}=\frac{1}{|\Lambda|}\sum_{\lambda\in\Lambda}F_{XY,\lambda}^{N}$. Consider the following two conditions for $\gamma<2$ and $\beta>0$:

[TABLE]

[TABLE]

[TABLE]

Then

[TABLE]

Lemma 18 is nearly identical to Lemma 7 in Krishnamurthy et al. [7] which is itself a modification of Theorem 2.2 in Tsybakov [125]. The main change is the addition of the second option to condition (ii). Thus the only required addition to the proof is to check that the second option is indeed sufficient. Since we consider the probability of the absolute value of $\hat{G}_{N}-I\left(f_{XY},f_{X}f_{Y}\right)$ being greater than $\beta$, it is clear that following the original proof with the second option will give the same result as the first option in condition (ii).

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

Lemma 19**.**

Let $R_{1},\dots,R_{\ell}$ be a partition of $[0,1]^{d_{X}+d_{Y}}$, each being cubes with side length $\ell^{-1/\left(d_{X}+d_{Y}\right)}$. Then there exists functions $u_{1},\dots,u_{\ell}$ such that,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $B(z,\epsilon)$ denotes an $L_{2}$ ball around z=\left[\begin{array}[]{c}x\\ y\end{array}\right] with radius \epsilon$$\in(0,1).

Proof:

Let $\epsilon>0$. We can construct two orthonormal systems of $q>3$ functions. Construct the first system on $[0,1]^{d_{X}}$ such that $\phi_{X,1}=1$, $\text{supp}\left(\phi_{X,j}\right)\subset[\epsilon,1-\epsilon]^{d_{X}}$, and $\left\|D^{r}\phi_{X,j}\right\|_{\infty}\leq J_{X}<\infty$ for all $j$. The second system is constructed on $[0,1]^{d_{Y}}$such that $\phi_{Y,1}=1$, $\text{supp}\left(\phi_{Y,j}\right)\subset[\epsilon,1-\epsilon]^{d_{Y}}$, and $\left\|D^{r}\phi_{Y,j}\right\|_{\infty}\leq J_{Y}<\infty$ for all $j$. We can then construct a combined orthonormal system $\phi_{i,j}=\phi_{X,i}\phi_{Y,j}$ which has $q^{2}$ functions. It is clear that this is an orthonormal system since $\int\phi_{i,j}(x,y)\phi_{m,n}(x,y)dxdy=1$ when $i=m$ and $j=n$ and zero otherwise. Also, $\phi_{1,1}=1$, $\text{supp}\left(\phi_{i,j}\right)\subset[\epsilon,1-\epsilon]^{d_{X}+d_{Y}},$ and there exists some $J<\infty$ such that $\left\|D^{r}\phi_{i,j}\right\|_{\infty}\leq J$ for all $i,j$.

Now for any given function $f\in L_{2}\left([0,1]^{d_{X}+d_{Y}}\right)$, we can find a unit-normed function $v\in\text{span}\left(\left\{\phi_{i,j}\right\}\right)$ such that $v\perp\phi_{1,i}$ for all $i$, $v\perp\phi_{i,1}$ for all $i$, and $v\perp f$. We can write $v=\sum_{i,j=1}^{q}a_{i}b_{j}\phi_{i,j}$. Then $D^{r}v=\sum_{i,j=1}^{q}a_{i}b_{j}D^{r}\phi_{i,j}$ which implies that

[TABLE]

where the last inequality comes from the fact that $v$ is unit-normed. Now define $\nu=\frac{1}{Jq}v$. Then clearly $\int\nu^{2}(x,y)dxdy$ is upper and lower bounded and we have that $\left\|D^{r}\nu\right\|_{\infty}\leq 1$.

We can now construct the functions $u_{j}$. First map $R_{j}$ to $[0,1]^{d_{X}+d_{Y}}$by scaling it appropriately. Then set

[TABLE]

where $\mathbf{j}$ is the point in $R_{j}$ that is mapped to $\mathbf{0}$ after scaling. This maps $[0,1]^{d_{X}+d_{Y}}$ back to $R_{j}$ while inheriting the properties derived from the construction of $\nu$. Now let $f$ be $g^{\prime}\left(\frac{f_{X}(x)f_{Y}(y)}{f_{XY}(x,y)}\right)\frac{f_{X}(x)f_{Y}(y)}{f_{XY}(x,y)}$ constrained to $R_{j}$ and scaled to fit $[0,1]^{d_{X}+d_{Y}}$. Conditions 1, 3, and 4 above are then fulfilled by construction. Also $\int_{R_{j}}u_{j}^{2}(x,y)dxdy=\frac{1}{\ell}\nu^{2}(x,y)dxdy\in\Theta\left(\ell^{-1}\right)$ which is condition 2 above. It’s also clear that $\left\|D^{r}u_{j}\right\|_{\infty}\leq\ell^{r/\left(d_{X}+d_{Y}\right)}$ which is condition 5, completing the proof. ∎

We now can prove Theorem 4. We will construct the conditions necessary to apply Lemma 18. Apply Lemma 19 to obtain an index set $\tilde{\Lambda}=\{-1,1\}^{\ell}$ and functions $u_{1},\dots,u_{\ell}$. Define the following set of perturbed functions around $f_{XY}$:

[TABLE]

This will form our set of alternatives. Due to the fact that $\int u_{j}(x,y)dx=\int u_{j}(x,y)dy=0$, we have that

[TABLE]

[TABLE]

That is, the perturbations on $f_{XY}$ are chosen so that the resulting marginal distributions are unperturbed.

The perturbation functions $u_{j}$ in Lemma 19 are restricted to the small $R_{j}$ bins and thus violate the Hölder class assumption. However, by scaling $H_{1}$ appropriately, we can ensure that $f_{XY,\lambda}\in\Sigma(s,H)$. We show this by following the same argument as in Krishnamurthy et al. [7], which we repeat here for completeness. Define $u_{\lambda}=H_{1}\sum_{j=1}^{\ell}\lambda_{j}u_{j}$. We will first show that $u_{\lambda}$ is Hölder smooth. Then $f_{XY,\lambda}$ is Hölder smooth by the triangle inequality. For $u_{\lambda}$, fix two points $v,z\in\mathbb{R}^{d_{X}+d_{Y}}$ and fix $r$ with $\sum_{j}r_{j}=s$. Define $z_{1}$ as the boundary point of the $R_{j}$ bin containing $z$ along the line between $z$ and $v$ and define $v_{1}$ similarly as the boundary point for the bin containing $v$ along the same line. We then have the following:

[TABLE]

The first line is an application of the triangle inequality. The second line follows from the fact that $u_{\lambda}$ and all of its derivatives are zero on the boundaries of the cubes $R_{j}$ as $u_{j}$ is not supported in the band around the border of $R_{j}$. The third line follows from the fundamental theorem of calculus where $\gamma(z,z_{1})$ is the path between $z$ and $z_{1}$. The fourth line is an application of Hölder’s inequality where we replace each derivative with its supremum, leaving just the path integral which simplifies to the length of the path. The fifth line follows from the assumption that $\left\|D^{r}u_{j}\right\|_{\infty}\leq\ell^{r/\left(d_{X}+d_{Y}\right)}$ when $\sum_{j}r_{j}\leq s+1$. For the sixth line, since $z$ and $z_{1}$ are in the same bin, then $\left\|z-z_{1}\right\|_{2}\leq\sqrt{d_{X}+d_{Y}}\ell^{-1/\left(d_{X}+d_{Y}\right)}$ as there are $\ell$ boxes with side length $\ell^{-1/\left(d_{X}+d_{Y}\right)}$. Finally, the last line follows since $z_{1}$ and $v_{1}$ are on the line segment between $z$ and $v$.

This indicates that $u_{\lambda}$, and therefore $f_{XY,\lambda}$, is guaranteed to be Hölder smooth if $H_{1}\ell^{\frac{s}{d_{X}+d_{Y}}}\sqrt{d_{X}+d_{Y}}\leq H$. Thus we require that $H_{1}=O\left(\ell^{-\frac{s}{d_{X}+d_{Y}}}\right)$. We will set $\ell$ later on.

Note that for any $f_{XY,\lambda}\in\Lambda$, by a second order Taylor series approximation in the first argument we have

[TABLE]

where $f_{XY}^{*}$ is the function from Taylor’s remainder theorem. By construction (see Lemma 19), the first order term vanishes. For the second order term, note that $f_{XY}^{*}$ lies on the line segment between $f_{XY}$ and $f_{XY,\lambda}$ and is therefore upper and lower bounded. Similarly the density $f_{XY,\lambda}$ will be upper and lower bounded for $N$ sufficiently large as $f_{XY,\lambda}\in\left[f_{XY}-H_{1},f_{XY}+H_{1}\right]$ due to the fact that $\left\|D_{0}u_{j}\right\|_{\infty}=\left\|u_{j}\right\|_{\infty}\leq 1,$ and $H_{1}$ will be chosen to decrease as $N$ increases. Assume without loss of generality that given $\epsilon>0$, $g^{\prime\prime}(\epsilon)>0.$ Thus there exists a constant $c_{0}$ such that

[TABLE]

where we have used the facts that the $u_{j}$ functions are orthogonal to each other and $\left\|u_{j}\right\|_{2}^{2}\in\Theta\left(\ell^{-1}\right)$. Therefore,

[TABLE]

providing us with the necessary separation of $2\beta$ where $\beta=c_{1}H_{1}^{2}/2$. Note that if $g^{\prime\prime}(\epsilon)<0$, we can simply consider $I\left(f_{XY},f_{X}f_{Y}\right)-I\left(f_{XY,\lambda},f_{X}f_{Y}\right)$ instead.

We now focus on bounding the squared Hellinger distance $H^{2}\left(\bar{F}_{XY}^{N}\times\left(F_{X}F_{Y}\right)^{N},F_{XY}^{N}\times\left(F_{X}F_{Y}\right)^{N}\right)$ where $\bar{F}_{XY}^{N}=\frac{1}{|\Lambda|}\sum_{\lambda\in\Lambda}F_{XY,\lambda}^{N}$. The Hellinger distance decomposes across product measures resulting in:

[TABLE]

To bound this, we will use the following result from Birge and Massart [74]:

Lemma 20**.**

Consider a set of densities $p$ and $p_{\lambda}=p\left(1+\sum_{j}\lambda v_{j}\right)$ for $\lambda\in\Lambda=\{-1,1\}^{\ell}$. Suppose that (i) $\left\|v_{j}\right\|_{\infty}\leq 1$, (ii) $\left\|1_{\{R_{j}^{C}\}}v_{j}\right\|_{1}=0,$ (iii) $\int pv_{j}=0,$ and (iv) $\int pv_{j}^{2}=\alpha_{j}>0$ all hold. Define $\bar{P}^{N}=\frac{1}{|\Lambda|}\sum_{\lambda\in\Lambda}P_{\lambda}^{N}$. Then

[TABLE]

To apply Lemma 20, define $v_{j}(z)=H_{1}u_{j}(z)/f_{XY}(z)$. Then $f_{XY,\lambda}=f_{XY}\left(1+\sum_{j}\lambda_{j}v_{j}\right)$. Requirements (i)-(iii) are immediately satisfied based on the properties of $u_{j}$ (see Lemma 19). Furthermore,

[TABLE]

for some constant $C$. Therefore

[TABLE]

Set $\ell=N^{\frac{2\left(d_{X}+d_{Y}\right)}{4s+d_{X}+d_{Y}}}$ resulting in $H_{1}=N^{-\frac{2s}{4s+d_{X}+d_{Y}}}.$ Then the Hellinger distance is bounded by a constant. Additionally, the error is larger than $\beta\in\Theta\left(N^{-\frac{4s}{4s+d_{X}+d_{Y}}}\right)$ allowing us to apply Lemma 18 when $s<\left(d_{X}+d_{Y}\right)/4$. Markov’s inequality then finishes the proof.

For $s>\left(d_{X}+d_{Y}\right)/4,$we get a lower bound of $O\left(N^{-1}\right)$ which is the parametric rate. In general, we cannot do any better than this [126, 7] thus establishing the lower bound in this regime. In particular, Krishnamurthy et al. [7] use a contradiction approach to establish this for divergence estimation which can be extended to the MI estimation problem.

Appendix G Theory for Mixed Random Variables

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

G-A Proof of Lemma 5

For (24), note that$\mathbf{N}_{xy}$ is a binomial random variable with parameter $f_{X_{D}Y_{D}}\left(x,y\right)$, $N$ trials, and mean $Nf_{X_{D}Y_{D}}\left(x,y\right)$. Thus (24) is the (potentially) fractional moment of a binomial random variable. By the generalized binomial theorem, we have that

[TABLE]

From [81], the $i$-th central moment of $\mathbf{N}_{xy}$ has the form of

[TABLE]

Combining this with (79) gives

[TABLE]

For (25), we apply a Taylor series expansion to obtain

[TABLE]

where we set $p=f_{X_{D}Y_{D}}\left(x,y\right)$, $p_{x}=f_{X_{D}}(x)$, and $p_{y}=f_{Y_{D}}(y)$ for notational convenience. By taking the expected value with respect to $\mathbf{N}_{x}$, $\mathbf{N}_{y}$, and $\mathbf{N}_{xy}$, we obtain

[TABLE]

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

G-B Proof of Theorem 6 (Bias)

For notational ease, let

[TABLE]

We have that

[TABLE]

We consider the second term in (81) first. A Taylor series expansion of $g\left(\mathcal{T}(\mathbf{X},\mathbf{Y})\times\frac{\mathbf{N}_{x}\mathbf{N}_{y}}{N\mathbf{N}_{xy}}\right)$ evaluated at $\mathcal{T}(\mathbf{X},\mathbf{Y})\times\frac{f_{X_{D}}\left(x\right)f_{Y_{D}}\left(y\right)}{f_{X_{D}Y_{D}}\left(x,y\right)}$ gives terms of the form of

[TABLE]

[TABLE]

where $i$ is a positive integer. For notational ease, set $p=f_{X_{D}Y_{D}}\left(x,y\right)$. By applying the binomial theorem and (24), we obtain

[TABLE]

Using a similar approach with (25), it can be shown that

[TABLE]

Thus the second term in (81) reduces to $O(1/N)$.

By conditioning on $\mathbf{X}_{1,D},\dots,\mathbf{X}_{N,D},\mathbf{Y}_{1,D},\dots,\mathbf{Y}_{N,D}$, the first term in (81) can be written as

[TABLE]

The conditional bias of $\tilde{\mathbf{G}}_{h_{X_{C}|x},h_{Y_{C}|y}}$ given $\mathbf{X}_{1,D},\dots,\mathbf{X}_{N,D},\mathbf{Y}_{1,D},\dots,\mathbf{Y}_{N,D}$ can be obtained from Theorem 2 as

[TABLE]

This expression provides the motivation for our choice of $\mathbf{h}_{X_{C}|x}$ and $\mathbf{h}_{Y_{C}|y}$. Since $\mathbf{h}_{X_{C}|x}\propto\mathbf{N}_{x}^{-\beta}$ and $\mathbf{h}_{Y_{C}|y}\propto\mathbf{N}_{y}^{-\alpha}$, then (23) gives terms with the form of $\mathbf{N}_{xy}\mathbf{N}_{x}^{-\beta i}\mathbf{N}_{y}^{-\alpha j}/N$ with $i+j\geq 1$. From Lemma 5, taking the expected value of these terms gives

[TABLE]

Similarly, taking the expectation of $\mathbf{N}_{xy}\mathbf{N}_{x}^{\beta d_{X}}\mathbf{N}_{y}^{\alpha d_{Y}}/N^{2}$ gives $O\left(N^{\beta d_{X}+\alpha d_{Y}-1}\right)$. Note that the polynomial terms of $\mathbf{N}_{x}\mathbf{N}_{y}/N^{2}$ and $\mathbf{N}_{xy}/N$ in the constants in (84) do not contribute to the bias rate as the $\mathbf{N}_{x}\mathbf{N}_{y}$ and $\mathbf{N}_{xy}$ terms in the numerator are cancelled by the $N^{2}$ and $N$ terms in the denominator, respectively, after taking the expectation. Combining all of these results completes the proof.

G-C Proof of Theorem 7 (Variance)

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

[TABLE]

Given all of the $\mathbf{X}_{i,D}$ and $\mathbf{Y}_{i,D}$ random variables, the estimators $\tilde{\mathbf{G}}_{h_{X_{C}|x},h_{Y_{C}|y}}$ are all conditionally independent since they use different sets of $\mathbf{X}_{i,C}$’s and $\mathbf{Y}_{i,C}$’s for each pair $(x,y)$. Thus from Theorem 3, we get

[TABLE]

Taking the expectation yields $O(1/N)$.

For the second term in (85), we know from (84) that

[TABLE]

Let $\mathbf{N}_{xy}^{\prime}$, $\mathbf{N}_{x}^{\prime}$, and $\mathbf{N}_{y}^{\prime}$ be independent and identically distributed realizations of $\mathbf{N}_{xy}$, $\mathbf{N}_{x}$, and $\mathbf{N}_{y},$ respectively. Then by the Efron-Stein inequality,

[TABLE]

where since $\mathbf{N}_{x}$, $\mathbf{N}_{y}$, and $\mathbf{N}_{xy}$ are not independent, we consider the effect of resampling all three simultaneously. Note that

[TABLE]

By Jensen’s inequality, we can consider separately each of the squared differences in (91). Then since $\left(\mathbf{N}_{xy},\mathbf{N}_{x},\mathbf{N}_{y}\right)$ is independent of $\left(\mathbf{N}_{xy}^{\prime},\mathbf{N}_{x}^{\prime},\mathbf{N}_{y}^{\prime}\right)$ and they are identically distributed, then the expected squared difference is proportional to the variance. For example, applying Lemma 5 gives

[TABLE]

By a similar procedure, we obtain

[TABLE]

Combining these results with (86) and (85) completes the proof.

Appendix H Proof of Theorem 9 (CLT)

We will first find the asymptotic distribution of

[TABLE]

By the standard central limit theorem [127], the second term converges in distribution to a Gaussian random variable with variance

[TABLE]

All that remains is to show that the first term converges in probability to zero as Slutsky’s theorem [128] can then be applied. Denote this first term as $\mathbf{W}_{N}$ and note that $\mathbb{E}\left[\mathbf{W}_{N}\right]=0$.

We will use Chebyshev’s inequality combined with the Efron-Stein inequality to bound the variance of $\mathbf{W}_{N}$. Consider the samples $\left\{\mathbf{Z}_{1},\dots,\mathbf{Z}_{N}\right\}$ and $\left\{\mathbf{Z}_{1}^{{}^{\prime}},\mathbf{Z}_{2},\dots,\mathbf{Z}_{N}\right\}$ and the respective sequences $\mathbf{W}_{N}$ and $\mathbf{W}_{N}^{{}^{\prime}}$. This gives

[TABLE]

Note that

[TABLE]

We will use the Efron-Stein inequality to bound $\mathbb{V}_{\mathbf{Z}_{1}}\left[g\left(\frac{\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X}_{1})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y}_{1})\nu_{1}\nu_{2}}{\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{X}_{1},\mathbf{Y}_{1})\nu_{3}}\right)\right]$. We thus need to bound the conditional expectation of the term

[TABLE]

where $\mathbf{Z}_{i}$ is replaced with $\mathbf{Z}_{i}^{{}^{\prime}}$ in the KDEs for some $i\neq 1$. Using similar steps as in Appendix E, we have that

[TABLE]

Then by the Efron-Stein inequality, $\mathbb{V}_{\mathbf{Z}_{1}}\left[g\left(\frac{\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X}_{1})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y}_{1})\nu_{1}\nu_{2}}{\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{X}_{1},\mathbf{Y}_{1})\nu_{3}}\right)\right]=O\left(\frac{1}{N}\right)$. Therefore

[TABLE]

A similar result holds for the $g\left(\frac{\tilde{\mathbf{f}}_{X,h_{X}}(\mathbf{X}_{1}^{{}^{\prime}})\tilde{\mathbf{f}}_{Y,h_{Y}}(\mathbf{Y}_{1}^{{}^{\prime}})\nu_{1}\nu_{2}}{\tilde{\mathbf{f}}_{Z,h_{Z}}(\mathbf{X}_{1}^{{}^{\prime}},\mathbf{Y}_{1}^{{}^{\prime}})\nu_{3}}\right)$ term in (92).

For the third term in (92),

[TABLE]

For the $N-1$ terms where $i=j$, we know from Appendix E that

[TABLE]

Thus these terms contribute $O(1/N)$. For the $(N-1)^{2}-(N-1)$ terms where $i\neq j$, we can do multiple substitutions of the form $\mathbf{u}_{j}=\frac{\mathbf{X}_{j}-\mathbf{X}_{1}}{h_{X}}$ resulting in

[TABLE]

Since $h_{X}^{d_{X}}h_{Y}^{d_{Y}}=o(1)$,

[TABLE]

Combining all of these results with Jensen’s inequality gives

[TABLE]

Applying the Efron-Stein inequality gives that $\mathbb{V}\left[\mathbf{W}_{N}\right]=o(1)$. Then by ChebyShev’s inequality, $\mathbf{W}_{N}$ converges to zero in probability. This completes the proof for the plug-in estimator.

For the weighted ensemble estimator, we present a more general result where we have different parameters $l_{X}\in\mathcal{L}_{X}$ and $l_{Y}\in\mathcal{L}_{Y}$ for $\mathbf{h}_{X_{C}|x}$ and $\mathbf{h}_{Y_{C}|y}$, respectively. We can then write

[TABLE]

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

[TABLE]

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

[TABLE]

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 10 follows a similar procedure as the extension to the ensemble case.