Direct estimation of density functionals using a polynomial basis

TL;DR

This paper introduces a novel data-driven basis approach for directly estimating density functionals from samples, avoiding distribution fitting and high-dimensional integration, with applications to divergence measures and error bounds.

Contribution

It proposes a new complete basis for approximating density functionals directly from samples, eliminating the need for distribution assumptions or explicit density estimation.

Findings

Accurate estimation of divergence functions from samples.

Effective approximation of density functionals using the new basis.

Empirical bounds on Bayes error rate demonstrated.

Abstract

A number of fundamental quantities in statistical signal processing and information theory can be expressed as integral functions of two probability density functions. Such quantities are called density functionals as they map density functions onto the real line. For example, information divergence functions measure the dissimilarity between two probability density functions and are useful in a number of applications. Typically, estimating these quantities requires complete knowledge of the underlying distribution followed by multi-dimensional integration. Existing methods make parametric assumptions about the data distribution or use non-parametric density estimation followed by high-dimensional integration. In this paper, we propose a new alternative. We introduce the concept of "data-driven basis functions" - functions of distributions whose value we can estimate given only samples from the underlying distributions without requiring distribution fitting or direct integration. We derive a new data-driven complete basis that is similar to the deterministic Bernstein polynomial basis and develop two methods for performing basis expansions of functionals of two distributions. We also show that the new basis set allows us to approximate functions of distributions as closely as desired. Finally, we evaluate the methodology by developing data driven estimators for the Kullback-Leibler divergences and the Hellinger distance and by constructing empirical estimates of tight bounds on the Bayes error rate.

Figures

Tables

Table 1. TABLE I: Experiment overview table.

	$f_0\left(𝐱\right)$			$f_1\left(𝐱\right)$
	Family	$\mu$	$\beta$	Family	$\mu$	$\beta$
Experiment 1	Normal	0	0	Normal	$\sqrt{\frac{1}{3}}$	0
Experiment 2	Normal	0	0.8	Normal	$\sqrt{\frac{1}{3}}$	0.8
Experiment 3	Normal	0	0.8	Normal	$\sqrt{\frac{1}{3}}$	0.9
Experiment 4	Normal	0	0	Uniform	0	3

Table 2. TABLE II: Parameters for 2 8-dimensional Gaussian data sets for which the Bayes error rate is known (from [ 57 ] ).

${𝒟}_0$	${𝝁}_{\mathrm{𝟎}}$	0	0	0	0	0	0	0	0
	${𝝈}_{\mathrm{𝟎}}$	1	1	1	1	1	1	1	1
	${𝝁}_{\mathrm{𝟏}}$	2.56	0	0	0	0	0	0	0
	${𝝈}_{\mathrm{𝟏}}$	1	1	1	1	1	1	1	1
${𝒟}_2$	${𝝁}_{\mathrm{𝟎}}$	0	0	0	0	0	0	0	0
	${𝝈}_{\mathrm{𝟎}}$	1	1	1	1	1	1	1	1
	${𝝁}_{\mathrm{𝟏}}$	3.86	3.10	0.84	0.84	1.64	1.08	0.26	0.01
	${𝝈}_{\mathrm{𝟏}}$	8.41	12.06	0.12	0.22	1.49	1.77	0.35	2.73

Table 3. TABLE III: Comparing upper bounds on the Bayes error rate for the multivariate Gaussians defined in Table II .

	Data 1	Data 2
Actual Bayes Error	10%	1.90%
Mahalanobis Bound	18.90% $\pm$ 0.55%	14.07 % $\pm$ 0.45 %
Bhattacharyya Bound	21.74% $\pm$ 0.87 %	4.68 % $\pm$ 0.27 %
$D_p$ Bound	16.51 % $\pm$ 1.07%	3.99 % $\pm$ 0.52 %
Convex Bound	14.17 % $\pm$ 0.86%	3.87% $\pm$ 0.43%

Equations

$$f_{\mathbf{x}}(\mathbf{x})=p_0{f}_0(\mathbf{x})+p_1{f}_1(\mathbf{x})$$

$$\eta ({\mathbf{x}}^{\ast })$$

$$=\frac{p_1{f}_1({\mathbf{x}}^{\ast })}{f_{\mathbf{x}}({\mathbf{x}}^{\ast })}=\frac{p_1{f}_1({\mathbf{x}}^{\ast })}{p_0{f}_0({\mathbf{x}}^{\ast })+p_1{f}_1({\mathbf{x}}^{\ast })}$$

$$P[y=0|\mathbf{x}={\mathbf{x}}^{\ast }]=\frac{p_0{f}_0({\mathbf{x}}^{\ast })}{p_0{f}_0({\mathbf{x}}^{\ast })+p_1{f}_1({\mathbf{x}}^{\ast })},$$

$$P[y=0|\mathbf{x}={\mathbf{x}}^{\ast }]=1-\eta ({\mathbf{x}}^{\ast }).$$

$$G(f_0,f_1)=\int g(\eta (\mathbf{x}))f_{\mathbf{x}}(\mathbf{x})d\mathbf{x}.$$

$$D_{\phi}(f_{0},f_{1})=\int\phi\bigg{(}\frac{f_{0}(\mathbf{x})}{f_{1}(\mathbf{x})}\bigg{)}f_{1}(\mathbf{x})d\mathbf{x},$$

$$\frac{f_0(\mathbf{x})}{f_1(\mathbf{x})}=\frac{p_1(1-\eta (\mathbf{x}))}{p_0\eta (\mathbf{x})}$$

$$f_1(\mathbf{x})=\frac{\eta (\mathbf{x})}{p_1}{f}_{\mathbf{x}}(\mathbf{x})$$

$$D_{\phi}(f_{0},f_{1})=\int\phi\bigg{(}\frac{p_{1}(1-\eta(\mathbf{x}))}{p_{0}\eta(\mathbf{x})}\bigg{)}\frac{\eta(\mathbf{x})}{p_{1}}f_{\mathbf{x}}(\mathbf{x})d\mathbf{x}.$$

$$g(\eta)=\frac{\eta}{p_{1}}\phi\bigg{(}\frac{p_{1}(1-\eta)}{p_{0}\eta}\bigg{)}.$$

$$H_i(f_0,f_1)=\int h_i(\eta (\mathbf{x}))f_{\mathbf{x}}(\mathbf{x})d\mathbf{x}.$$

$$g(\eta )\approx \sum _{i=0}^k{w}_i{h}_i(\eta ),$$

$$G(f_0,f_1)\approx \hat{G}(f_0,f_1)=\sum _{i=0}^k{w}_i{H}_i(\eta ),$$

$${\Phi }_k({\mathbf{x}}^{\ast })=\sum _{i:{\mathbf{x}}_i\in {\mathcal{N}}_k({\mathbf{x}}^{\ast })}{y}_i.$$

$$I_{r,k}(\mathbf{x})=\{\begin{array}{ll}1& {\Phi }_k(\mathbf{x})=r\\ 0& otherwise,\end{array}$$

$${\rho }_{r,k,N}(\mathbf{X})=\frac{1}{N}\sum _{\mathbf{x}\in \mathbf{X}}{I}_{r,k}(\mathbf{x}).$$

$$\underset{N\to \infty }{\lim }{\rho }_{r,k,N}(\mathbf{X})\stackrel{L^2}{\to }\int \left(\genfrac{}{}{0ex}{}{k}{r}\right){\eta }^r(\mathbf{x})(1-\eta (\mathbf{x}){)}^{k-r}{f}_{\mathbf{x}}(\mathbf{x})d\mathbf{x}$$

$$H_{r,k}(f_0,f_1)$$

$$h_{r,k}(\eta )=\left(\genfrac{}{}{0ex}{}{k}{r}\right){\eta }^r(1-\eta {)}^{k-r}.$$

$$g(\eta)\approx\sum_{r=0}^{k}g\Big{(}\frac{r}{k}\Big{)}h_{r,k}(\eta).$$

$$G(f_0,f_1)=\int g(\eta (\mathbf{x}))f_{\mathbf{x}}(\mathbf{x})d\mathbf{x},$$

$$\hat{G}_{k,N}(\mathbf{X})=\sum_{r=0}^{k}g\Big{(}\frac{r}{k}\Big{)}\rho_{r,k,N}(\mathbf{X})$$

$$\lim\limits_{k\rightarrow\infty}\lim_{\begin{subarray}{c}N\to\infty\\ k/N\to 0\end{subarray}}E\bigg{[}\Big{(}\hat{G}_{k,N}(\mathbf{X})-G(f_{0},f_{1})\Big{)}^{2}\bigg{]}=0.$$

$${\hat{G}}_{k,N}(\mathbf{X})=\sum _{r=0}^k{w}_r{\hat{H}}_{r,k,N}(\mathbf{X}).$$

$$e_T$$

$$=\underset{=:e_A}{\underbrace{G(f_0,f_1)-{\hat{G}}_k(f_0,f_1)}}+\underset{=:e_{est}}{\underbrace{{\hat{G}}_k(f_0,f_1)-{\hat{G}}_{k,N}(\mathbf{X})}}.$$

$$g(\eta )=\left(\genfrac{}{}{0ex}{}{3}{1}\right)\eta (1-\eta {)}^2$$

$$\sum _{r=0}^3{w}_r{\beta }_{r,3}(\eta )=g(\eta ),$$

$$\hat{g}(\eta )=$$

$$=$$

$$\displaystyle+g\Big{(}\frac{3}{3}\Big{)}\beta_{3,3}(\eta)$$

Full text

I Introduction

Information divergence measures play a central role in the fields of machine learning and information theory. Information divergence functions, functionals that map density functions to IR, have been used in many signal processing applications involving classification [1], segmentation [2], source separation [3], clustering [4], and other domains. In machine learning, a sub-class of these divergences known as $f$-divergences [5], are widely used as surrogate loss functions since they form convex upper bounds on the non-convex 0-1 loss [6].

Although these measures prove useful in a variety of applications, the task of estimating them from multivariate probability distributions using finite sample data can pose a significant challenge. It is especially challenging for continuous distributions, which is the focus of this paper. For algorithms and theory for estimating information measures for discrete distributions the reader is referred to [7, 8, 9, 10, 11, 12]. For the continuous case treated here, there are three general classes of methods for estimating divergence [13]: 1) parametric methods, 2) non-parametric methods based on density estimation, and 3) non-parametric methods based on direct (or graph-based) estimation. Parametric methods are the most common choice for estimation, and typically offer good convergence rates ($1/N$) when an accurate parametric model is selected. The fundamental limitation of parametric methods, is that an accurate parametric model is rarely available in real world problems and using an inaccurate parametric model can heavily bias the final estimate. As an alternative, when no parametric form is known, non-parametric density estimates such as kernel density estimation [14], histogram estimation [15], or $k$-nearest neighbor ($k$-NN) density estimation [16] are used to characterize the distribution. While these methods are quite powerful in certain scenarios, they are generally high variance, sensitive to outliers, and scale poorly with dimension [13].

An alternative to these two classes of methods is direct (or graph-based) estimation, which exploits the asymptotic properties of minimal graphs in order to directly estimate distribution functionals without ever estimating the underlying distributions themselves. These methods have been used to estimate density functionals such as entropy [17, 18, 19], the information divergence [13], and the $D_{p}$-divergence [20]. This class of methods can have faster asymptotic convergence rates [13] and are often simpler to implement than plug-in methods which may have many tuning parameters such as kernel width or histogram bin size. Direct estimators are often inspired by an asymptotic property of a graph-theoretic quantity that can be scaled or modified in order to generate an estimator for a given information-theoretic measure. This approach is customized to a specific form of the density functional and is sometimes difficult to generalize. The estimator for the $D_{p}$ divergence is an example of this approach [20]. This is in contrast with plug-in estimators, where the same general approach can be used to estimate any distribution functional. Estimators based on influence functions attempt to bridge the gap between the two types of approaches [21]. In that work the authors present a recipe for estimating any smooth functional by using a Von Mises expansion - the analog of the Taylor expansion for distributions; however that approach still requires that part of the data be used for density estimation.

In this paper, we provide a general approach for estimating a wide range of distribution functionals. We propose decomposing the functional onto a complete set of “data-driven” basis functions; where the term “data-driven” means that the basis is determined directly from the data and does not involve distribution fitting or direct integration. We show that a broad class of distribution functionals can be approximated as closely as desired through linear combinations of our proposed basis, where the weights of the basis expansion are determined through convex optimization. Our approach offers a powerful alternative for estimating information-theoretic distribution functionals. We demonstrate the flexibility of the approach by constructing empirical estimators of bounds on the Bayes error rate using the same basis set.

The remainder of the paper is organized as follows. In the next section, we review the literature in this area. In Section II we provide a detailed description of the problem this paper attempts to solve and establish some of the basic mathematical notation used throughout this paper. In Section III we introduce a set of graph-theoretic basis functions and prove that a wide range of information theoretic quantities can be represented by a linear combination of functions in this set. In sections IV and V, we explore the limitations of the proposed methodology in the finite sample regime and propose two alternate fitting routines to identify weights to map these basis functions to quantities of interest. In section VI we empirically investigate how the proposed method can be used to estimate popular divergence measures (the KL-divergence, the Hellinger distance, the $D_{p}$-divergence), and we compare its performance to various parametric and non-parametric alternatives. In Section VII, we show how the method can be extended to form tighter bounds on the Bayes error rate for binary classification problems. Section VIII offers some concluding remarks.

I-A Related Work

A natural method for non-parametric estimation of continuous distribution functionals involves histogram binning followed by plug-in estimation [22, 23]. When the bin-size is adjusted as a function of the number of available samples per bin, this histogram plug-in method is known as Grenander’s method of sieves and it enjoys attractive non-parametric convergence rates [24, 25]. While these methods may work well for small data dimension ($d=1,2$), their complexity becomes prohibitive for larger dimensions. Recent work has focused on non-parametric plug-in estimators that are more practical in higher dimensions. These approaches generally only estimate the PDF for the values of samples in the reference data, then calculate the expected value across the sample data in place of numerical integration [26, 27, 28]. To circumvent the slow convergence associated with these approaches, ensemble methods [27] and methods based on influence functions [21] have been proposed which are capable of achieving the parametric rate $\mathcal{O}(N^{-1})$ MSE convergence if the underlying densities meet certain smoothness conditions [27]. Alternatively, estimates of divergence functions that rely on estimates of the likelihood ratio instead of density estimation have been proposed for estimating the $\alpha$-divergence and the $L_{2}$-divergence [6, 29, 30, 31, 32]. These methods estimate the likelihood ratio of the two density functions and plug that value into the divergence functions. Other approaches that bypass density estimation are the convex optimization approach of [29] to estimate $f$-divergences and the $k$-NN graph and minimal spanning tree approaches to estimating Henze-Penrose divergence [33, 34, 20].

Similarly, the approach we propose in this paper bypasses density estimation through a polynomial basis expansion where the basis coefficients are determined through a convex optimization criterion. This provides added flexibility and allows us to easily estimate a large class of distribution functionals and to establish empirical estimates of bounds on the Bayes error rate. Conceptually, this approach is similar to prior work in estimating the entropy of discrete distributions using polynomial approximations [7, 8, 9, 11].

Bounds on optimal performance are a key component in the statistical signal processing literature. For classification problems, it is often desirable to bound the Bayes error rate (BER) - the minimum achievable error in classification problems. The well-known Chernoff upper bound on the probability of error has been used in a number of statistical signal processing applications [35]. It motivated the Chernoff $\alpha$-divergence [13]. The Bhattacharyya distance, a special case of the Chernoff $\alpha$-divergence for $\alpha=\frac{1}{2}$, upper and lower bounds the BER [36, 37]. Beyond the bounds on the BER based on divergence measures, a number of other bounds exist based on other functionals of the distributions [38, 39]. For estimation problems, the Fisher information matrix (FIM) bounds the variance of the optimal unbiased estimator (through it’s relationship with the CRLB). The authors have also previously introduced the $D_{p}$ divergence, a non-parametric $f$-divergence, and showed that it provides provably tighter bounds on the BER than the BC bound [20]. They extended this work to estimation of the Fisher information in [40].

Our data-driven basis, consisting of Bernstein polynomials, can be used to estimate functionals of distributions and to estimate bounds on Bayes optimal classification performance. Bernstein polynomials of a different form have been used for density estimation [41, 42, 43, 44, 45, 46]. In contrast to this work, our methods do not rely on density estimation.

II Problem Setup

In this section, we will set up the problem and establish the notation that will be used throughout the rest of the paper. We are given a set of data $[\mathbf{X},\mathbf{y}]$ containing $N$ instances, where each instance is represented by a $d$-dimensional feature vector $\mathbf{x}_{i}$ and a binary label $y_{i}$. Suppose that this data is sampled from underlying distribution, $f_{\mathbf{x}}(\mathbf{x})$, where

[TABLE]

is made up of the two conditional class distributions $f_{0}(\mathbf{x})$ and $f_{1}(\mathbf{x})$ for classes [math] and $1$, with prior probabilities $p_{0}$ and $p_{1}$ respectively. If the priors aren’t explicitly known, they can be easily estimated from the sample data by measuring the ratio of samples drawn from each class. As a simple application of Bayes theorem, we can define the posterior likelihood of class 1, $\eta(\mathbf{x})$, evaluated at a point $\mathbf{x}=\mathbf{x}^{*}$, as

[TABLE]

We can similarly define the posterior probability for class 0 as

[TABLE]

and since $y$ is binary,

[TABLE]

To simplify the notation, we remove the dependence of $\eta$ on $\mathbf{x}^{*}$ from portions of the analysis that follow.

Suppose that we wish to estimate some functional $G(f_{0},f_{1})$ of distributions $f_{0}(\mathbf{x})$ and $f_{1}(\mathbf{x})$, which can be expressed in the following form

[TABLE]

Throughout the rest of this paper, we will refer to the $g(\eta)$ in (5) as the posterior mapping function. Many functionals in machine learning and information theory, such as $f$-divergences and loss functions, can be expressed this way. Consider the family of $f$-divergences as an example. They are defined as

[TABLE]

where $\phi(t)$ is a convex or concave function unique to the given $f$-divergence. By substituting

[TABLE]

and

[TABLE]

we can redefine (6) as

[TABLE]

Thus any $f$-divergence can be presented in the form outlined in (5) simply by defining the posterior mapping function $g(\eta)$ as

[TABLE]

We propose a procedure for estimating these types of divergence functionals which bypasses density estimation. We do this by representing the functional in terms of the asymptotic limit of a linear combination of graph-theoretic basis functions.

Suppose that there exists a set of basis functions $H_{0}(\eta),...,H_{k}(\eta)$ that can be similarly expressed as

[TABLE]

If we assume that there exists a set of coefficients such that

[TABLE]

then consequently

[TABLE]

where the sense of approximation is that the $\ell_{2}$ norm of the difference between the right and left hand sides is small. In the following section, we will introduce a set of basis functions that have the desired properties.

III Graph-theoretic Basis Functions

Consider the dataset $[\mathbf{X},\mathbf{y}]$ previously defined. Suppose we select an arbitrary instance $\mathbf{x}^{*}$ from $\mathbf{X}$ and examine it along with the set of its $k-1$ nearest neighbors $\mathbf{x}_{NN}^{1},\mathbf{x}_{NN}^{2},...,\mathbf{x}_{NN}^{k-1}$ in $\mathbf{X}$. We can define the neighborhood set $\mathcal{N}_{k}(\mathbf{x}^{*})=[\mathbf{x}^{*},\mathbf{x}_{NN}^{1},...,\mathbf{x}_{NN}^{k-1}]$, as the union of $\mathbf{x}^{*}$ and its $k-1$ nearest neighbors. Using this we define $\Phi_{k}(\mathbf{x}^{*})$ as the number of instances in the neighborhood set which are drawn from class 1, or alternatively, the sum of $y$ across all points in the neighborhood

[TABLE]

Figure 1 provides a simple illustration to help explain how $\Phi_{k}(\mathbf{x}^{*})$ is calculated. Calculating $\Phi_{k}$ is similar to how nearest neighbor classifiers make decisions, but with two important differences:

The base instance $\mathbf{x}^{*}$ is considered in the neighborhood indistinguishably from other instances in $\mathcal{N}_{k}(\mathbf{x}^{*})$. 2. 2.

Where traditional $k$-NN classifiers are concerned only with identifying the majority, we are interested in the exact number of instances drawn from each class.

In essence $\Phi_{k}(\mathbf{x}^{*})$ tells us something about the probability that $y=1$ for instances on or near $\mathbf{x}^{*}$. Since we are more concerned with the dataset as a whole than the local characteristics in $\mathbf{x}$, we define the statistic $\rho_{r,k,N}(\mathbf{X})$ to be the fraction of instances $\mathbf{x}\in\mathbf{X}$, for which $\Phi_{k}(\mathbf{x})=r$, $r\leq k$. If we define the indicator function $I_{r,k}(\mathbf{x})$ as

[TABLE]

then this test statistic $\rho_{r,k,N}(\mathbf{X})$ can be represented by

[TABLE]

The function $\rho_{r,k,N}(\mathbf{X})$ is simply the proportion of $k$-NN neighborhoods that contain exactly $r$ points from class $y=1$. This statistic has a number of desirable qualities. We show that this statistic asymptotically converges to a function of the underlying distributions that can be described in the form outlined in (5). The following is proven in Appendix A.

Theorem 1**.**

As the number of samples ($N$) approaches infinity,

[TABLE]

whenever $k/N\rightarrow 0$.

We propose to use the asymptotic form of $\rho_{r,k,N}(\mathbf{X})$ defined in Theorem 1 as a basis function for estimating functionals of the form (5),

[TABLE]

where

[TABLE]

The function (18) is the $r^{th}$ Bernstein basis polynomial of degree $k$ [47]. Bernstein’s proof of the Weierstrass Approximation Theorem [48] asserts that any continuous function $g(\eta)$ can be uniformly approximated on $\eta\in[0,1]$ to any desired accuracy by a linear combination of functions in (18) of the form

[TABLE]

Combining this result with Theorem 1, we can show that a linear combination of this basis can be used to estimate any function of the form (5) .

Theorem 2**.**

For any $G(f_{0},f_{1})$ that can be expressed in the form

[TABLE]

where $g(\eta)$ is continuous on $[0,1]$, the approximation

[TABLE]

satisfies

[TABLE]

Theorem 21 provides an asymptotically consistent method of estimating a variety of information-theoretic functions that makes no assumptions on the underlying distributions and can be calculated without having to perform density estimation. Throughout the rest of the paper we will refer to the weights $g(r/k)$ in the approximation specified by (20) in Theorem 21 as the Bernstein weights. We next turn to the finite sample properties of the estimator $\rho_{r,k,N}(\mathbf{X})$.

IV Finite Sample Considerations

The previous Section investigated the asymptotic properties of linear combinations of the proposed set of empirically estimable basis functions. The asymptotic consistency of the proposed method is valuable, however in real world scenarios, data is inherently a finite resource, and as a result the efficacy of this method is heavily dependent on its convergence characteristics in the finite sample regime. In this section, we will take a detailed look into how restricting both $N$ and $k$ affects our ability to estimate functions of two distributions. To do this, it is necessary to first break down the different sources of error in the proposed methodology.

IV-A Estimation vs. Approximation Error

The goal of this paper is to empirically estimate the functional $G(f_{0},f_{1})$ of the two underlying distributions $f_{0}(\mathbf{x})$ and $f_{1}(\mathbf{x})$ using a linear combination of directly estimable basis functions

[TABLE]

We divide the error of this estimate into two types, the approximation error ($e_{A}$) and the estimation error ($e_{est}$):

[TABLE]

This allows us to separate the error in estimating the basis functions from error in fitting to the posterior mapping function. Understanding the trade-off between these two error types will be particularly useful in Section V, where we explore different methods of fitting weights to the desired density functionals.

IV-B Considerations for finite $k$

A finite sample also implies a finite $k$ and impacts the approximation error. Let us consider the Bernstein weighting scheme introduced in (20) for the scenario where the size of the basis set ($k$) is restricted. Consider the following example problem.

Example: Suppose that we wish to estimate the function

[TABLE]

using the basis set $\beta_{0,3}(\eta),\beta_{1,3}(\eta),\beta_{2,3}(\eta),\beta_{3,3}(\eta)$. Because $g(\eta)=\beta_{1,3}(\eta)$, there exists a set of weights such that

[TABLE]

however, using the Bernstein weighting scheme in (20) yields

[TABLE]

It is clear from this example that the Bernstein weighting procedure do not always provide ideal weights when $k$ is restricted. Based on these results, we are motivated to explore alternative weighting procedures in order to improve the performance of this method for the finite sample case. In the following Section, we will introduce a method of finding better weights using convex optimization.

Regardless of how weights are assigned to the approximation, selection of $k$ remains an important factor affecting performance. In order to satisfy Theorem 21, $k$ should be functionally dependent on $N$, such that $N$ approaches infinity, both $k\rightarrow\infty$ and $\frac{k}{N}\rightarrow 0$, however this still provides a great degree of freedom in the selection of $k$.

In general, there are two major competing factors that must be considered when selecting $k$. The first is that the Weierstrass approximation theorem can exactly represent any posterior mapping function $g(\eta)$ as a linear combination of the proposed basis set only as $k\rightarrow\infty$. This provides motivation for choosing a large $k$-value to ensure the best possible fit of $g(\eta)$. The second factor is that the asymptotic characteristics of $\rho_{r,k,N}(\mathbf{X})$ are dependent on all points in $\mathcal{N}_{k}(\mathbf{x}^{*})$. Moreover the regime for which Theorem 21 holds, requires that $\frac{k}{N}\rightarrow 0$, so we are motivated to select $k$ such that $N\gg k$. This means that the selection of $k$ must achieve a compromise in the trade-off between the approximation and estimation errors, since larger values of $k$ will increase the amount of finite sample error made in estimating the individual basis functions, while lower values of $k$ may inhibit our ability to accurately model the desired function in the asymptotic regime.

To illustrate how our ability to estimate the desired set of basis functions varies with $k$, we calculate the estimated and asymptotic values of $\rho_{r,k,N}(\mathbf{X})$ for $k=10$ and $k=100$ on data drawn from distributions $f_{0}(\mathbf{x})\sim\mathcal{N}(\mathbf{0}_{3},\mathbf{I}_{3})$ and $f_{1}(\mathbf{x})\sim\mathcal{N}(\frac{1}{\sqrt{3}}\mathbf{1}_{3},\mathbf{I}_{3})$, where $\mathbf{0}_{3}=[0\ \ 0\ \ 0]$ and $\mathbf{1}_{3}=[1\ \ 1\ \ 1]$, and plot the results in Figure 2. Estimates are calculated at 3 different sample sizes ($N_{0}=N_{1}=100,1000,10000$) and each estimate shown in Figure 2 has been averaged across 500 Monte Carlo trials. While we can estimate the basis set for either $k$-value with a high degree of accuracy given enough samples, the estimates for $k=10$ are noticeably more accurate. In fact, we are able to do about as well with $1000$ samples for $k=10$ as we are with $10000$ samples for $k=100$.

V Optimization Criteria for Fitting Density Functionals

In this section, we propose a convex optimization criterion to identify appropriate weights for fitting information-theoretic functions when $k$ and $N$ are restricted. Inherently, our goal is to minimize the total error, defined in (23), however minimizing this quantity directly isn’t feasible since the value of $G(f_{0},f_{1})$ is unknown. To circumvent this challenge we focus on developing a criterion to minimize the approximation error. We initially develop an optimization criterion that assumes the posterior is uniformly distributed, then propose an alternate method which incorporates an estimate of the posterior density function in order to more accurately model the approximation error.

V-A Uniform Optimization Criteria

Recall that the approximation error $e_{A}$ can be represented as

[TABLE]

where

[TABLE]

Since solving (27) requires high-dimensional integration and knowledge of the underlying distributions. However, because $\eta$ is a function of $\mathbf{x}$, $\epsilon(\eta)$ is implicitly a function of $\mathbf{x}$ as well, and by the law of the unconscious statistician [49],

[TABLE]

where $f_{\eta}(\eta)$ is the probability density function of the random variable $\eta$. Rewriting the error in this form simplifies the region of integration to a well defined space (since $\eta\in[0,1]$) and circumvents the high dimensionality of $\mathbf{x}$. While this eliminates some of the challenges in calculating the error it also creates new ones stemming from the fact that $f_{\eta}(\eta)$ is unknown and difficult to estimate due to its implicit dependency on $f_{0}(\mathbf{x})$ and $f_{1}(\mathbf{x})$. The task of estimating $f_{\eta}(\eta)$ will be explored in detail in Section V-B, however for the time being we will bypass this challenge and simply attempt to minimize

[TABLE]

It is worth noting that if $f_{\eta}(\eta)$ is uniformly distributed

[TABLE]

While there exists an analytical solution for identifying the weights which minimize (30)[50], we use a convex optimization procedure that allows us to also account for the estimation error. If we define a discretized set of posterior values $\bm{\tilde{\eta}}=[\tilde{\eta}_{1},\tilde{\eta}_{2},...,\tilde{\eta}_{\tilde{N}}]$, where $0\leq\tilde{\eta}_{1}<\tilde{\eta}_{2}<...<\eta_{\tilde{N}}\leq 1$, a procedure to identify weights that minimize (30) can be defined as

[TABLE]

.

To illustrate the effectiveness of this method, we consider the example problem of trying to estimate the Hellinger distance (a problem we will further explore in Section VI). If we assume both classes have equal prior probability ($p_{0}=p_{1}=0.5$), then the posterior mapping function for the squared Hellinger distance is

[TABLE]

This function is estimated using this convex weighting procedure as well as the previously described Bernstein weighting procedure, and we compare how well each method models the desired function for values of $k$ varying from [math] to $100$. The performances of each method is evaluated by the following formula

[TABLE]

and the results are presented in Figure 3 for a range of $k$ values varying from $1$ to $100$. This experiment shows that the proposed convex fitting procedure is able to approximate the Hellinger posterior mapping function far more accurately than the Bernstein approximation. This improvement isn’t surprising since the proposed method directly minimizes the MSE whereas the Bernstein approximation guarantees consistency only as $k\rightarrow\infty$, but it still helps to illustrate the potential for improvement in the Bernstein weights that exists for smaller $k$.

The expression (34) does not take into account finite sample errors that lead to noisy estimates of the basis functions and thus does not directly reflect our ability to estimate $G(f_{0},f_{1})$ with a finite sample. The consequences of this could be quite significant. When using a similar approach for entropy estimation, Paninski observed that the good approximation properties were a result of “large oscillating coefficients” which magnify the variance of the corresponding estimator[11]. Additionally, it does not account for the possibility that $\eta$ is distributed non-uniformly.

To empirically examine the finite sample properties of the two approaches, we repeat the previous experiment, this time estimating the basis functions empirically on samples of data drawn from distributions $f_{0}(\mathbf{x})\sim\mathcal{N}(\mathbf{0}_{3},\mathbf{I}_{3})$ and $f_{1}(\mathbf{x})\sim\mathcal{N}(\frac{\mathbf{1}_{3}}{\sqrt{3}},\mathbf{I}_{3})$. We generate $N=1000$ samples (500 samples per class) in each of the 500 iterations of a Monte Carlo simulation, and evaluate the MSE as

[TABLE]

where $N_{MC}$ represents the number of Monte Carlo iterations. Since we know that the estimation error is scaled by the magnitude of the weights, we also evaluate a modified fitting routine which augments (32) with a regularization term to penalize solutions with large weights,

[TABLE]

where $\lambda$ represents a tuning parameter which controls the importance assigned to minimization of the approximation error relative to the estimation error. Intuitively, higher $\lambda$ values will make sense for smaller data sets to control the variance of the estimator. We set $\lambda=0.01$ for all experiments conducted in this paper. The results of this experiment are shown in Figure 4. We immediately see the necessity of the regularization term, as without it the error becomes extremely large for a range of $k$ values. More generally, the inclusion of the regularization term improves the performance at every $k$ value in this experiment. In comparing the Bernstein weights with the convex (regularized) weights, we find that 1) the performance of the convex method is less dependent on the selection of $k$ and 2) the convex weights generally perform better at lower values of $k$, while the Bernstein weights outperform at higher $k$. While the peak performance of the Bernstein method is higher than the convex method, there exists no good method of selecting $k$ a priori in order to reliably achieve this performance. In contrast, the convex method with regularization is less sensitive to the value of $k$ selected.

V-B Density-weighted Optimization Criteria

In the optimization criteria introduced in the previous section we implicitly make the assumption that the distribution of the random variable $\eta(\mathbf{x})\sim f_{\eta}(\eta)$ is uniformly distributed. In this section we will investigate a data-driven estimator for $f_{\eta}(\eta)$. However, before we proceed it is important to clarify what this distribution actually is.

We initially introduced $\eta$ as the posterior likelihood function for class 1, which we showed in (2) can be represented as a function of the underlying distributions. When this function’s input is a known point $\mathbf{x}^{*}$, $\eta(\mathbf{x}^{*})$ represents the probability that $y=1$ given that $\mathbf{x}=\mathbf{x}^{*}$. However, if the input is a random variable $\mathbf{x}$, then $\eta$ is also a random variable, which is distributed according to $f_{\eta}(\eta)$.

Figure 5 illustrates $f_{\eta}(\eta)$ for two univariate normal distributions $f_{0}(x)\sim N(0,1)$ and $f_{1}(x)\sim N(1,1)$. Figure 5a displays the two class distributions across $x$, Figure 5b displays $\eta(x)$ as a function of $x$, and Figure 5c displays $f_{\eta}(\eta)$ as a function of $\eta$. We see from this illustration that while, $\eta(x)$ is close to [math] or $1$ across most of the region of $x$ that is displayed, $\eta(x)$ remains close to $0.5$ in the regions where $f_{x}(x)$ is greatest. As a results $f_{\eta}(\eta)$ is roughly bell-shaped, and the likelihood of $\eta$ existing at the extremities (close to [math] or $1$) is relatively low. Because the probability density in these regions is low, the accuracy of our fit in these regions is less important, and can be given less weight in the fitting routine. Figure 6 repeats this illustration for two well seperated normal distributions. In this case the distribution of $\eta$ is such that $f_{\eta}(\eta)$ is most dense towards the extremities, and therefore they should be given more weight in the fitting routine. Side by side, these two figures present an interesting contrast. Despite the fact that the set of density functions $f_{0}$ and $f_{1}$ look quite similar in the two scenarios, the minor difference in separation significantly alters the distribution of $\eta$.

In practice the underlying distributions $f_{0}$ and $f_{1}$ are unknown, and as a result, $f_{\eta}(\eta)$ is unknown as well. However, if we were able to sample $f_{\eta}(\eta)$ at $\eta=\tilde{\eta}_{i}$, a more direct method of minimizing $e_{T}$ would be to solve

[TABLE]

where $\Delta_{\tilde{\eta}}=\tilde{\eta}_{i+1}-\tilde{\eta}_{i}$.

Below we show that $\rho_{r,k,N}(\mathbf{X})$, the statistic previously defined in Theorem 1 can be used to sample the PDF of $\eta$. This result is stated in Theorem 3.

Theorem 3**.**

As $N\rightarrow\infty$ and $k\rightarrow\infty$ in a linked manner such that $\frac{k}{N}\rightarrow 0$ and $\frac{r}{k}\rightarrow\eta^{*}$

[TABLE]

Theorem 3 is useful as it provides a method of sampling $f_{\eta}(\eta)$ that doesn’t depend on estimates of the underlying density functions $f_{0}$ and $f_{1}$. Using this result, we can estimate the density of the posterior at point $\tilde{\eta}_{i}$ as

[TABLE]

where $\tilde{r}_{i}$ and $\tilde{k}_{i}$ are integers selected such that $\tilde{\eta}_{i}=\frac{\tilde{r}_{i}}{\tilde{k}_{i}}$ for $i=1,2,...,\tilde{N}$. Sampling at exactly $\tilde{\eta}_{i}$ may not always be possible since the maximum value of $k$ is limited by the sample size, and $k$ determines the resolution of the sampling scheme. Even if it is possible, it may not be desirable to recalculate $\rho_{\tilde{r}_{i},\tilde{k}_{i},N}(\mathbf{X})$ for different values of $k_{i}$ because of the computational burden. To overcome these problems we can design our approach such that we utilize the same set of test statistics $\rho_{r,k,N}(\mathbf{X})$ in the estimation of the posterior distribution as are used in the estimation of the basis set. One way to do this is to assign the set of discretized posteriors

[TABLE]

so that it is straightforward to calculate $\hat{f}_{\eta}(\eta)$ from the known values of $\rho_{r,k,N}(\mathbf{X})$. An alternate approach is to leave $\tilde{\eta}$ unchanged and interpolate $\hat{f}_{\eta}(\eta)$ to determine its value at the desired points. The advantage of this approach is that it doesn’t constrain how $\eta$ is sampled. Throughout the rest of this paper, we will employ the latter method and solve for $\hat{f}_{\eta}(\eta)$ by linearly interpolating between its values on the discretized set (39).

Because this density-weighted fitting routine more directly minimizes the approximation error of the final estimate, we expect it to generally outperform the uniform method, particularly in cases where the density of the posterior is highly non-uniform and where the desired posterior mapping function $g(\eta)$ is difficult to model using the proposed basis set. This hypothesis will be verified in Sections VI and VII, when we empirically evaluate our methods with real data. This improvement comes at a computational cost since the weights are now data-dependent, they must be calculated online, whereas for the uniform method they can be calculated offline and stored. Solving for the $k$ weights can be be implemented with $\mathcal{O}(k^{3})$ time complexity [51]. Current state-of-the-art algorithms for $k$-NN graphs can be implemented with $\mathcal{O}(N\log N)$ time complexity [52]. As a result the complexity for the Convex (uniform) method is $\mathcal{O}(N\log N)$ assuming that the weights are available. The complexity for the Convex (density weighted) method is $\mathcal{O}(k(N)^{3}+N\log N)$ since the weights must be determined for each new dataset.

VI Divergence Estimation

In this section, we show how the proposed method can be applied to estimating three $f$-divergences, the Hellinger distance, the KL-divergence, and the $D_{p}$-divergence, directly from data.

VI-A Hellinger Distance

The Hellinger distance squared is an $f$-divergence used to quantify the dissimilarity between two probability distributions and is calculated by

[TABLE]

Using the approach proposed in Section V, we can estimate $H^{2}(f_{0},f_{1})$ by fitting weights to the posterior mapping function

[TABLE]

To evaluate the efficacy of this method, we conduct four different experiments in which we attempt to estimate the Hellinger distance between two distributions from finite sample data. In the first three experiments, both distributions are normally distributed according to $f(\mathbf{x})\sim N(\mu\mathbf{1}_{d},\mathbf{\Sigma}_{d})$, where

[TABLE]

for $\sigma_{i,j}=\beta^{|i-j|}$. The first experiment evaluates the most basic case where both Gaussians have spherical covariance. The second experiment considers the case where there exists a strong fixed dependency between adjacent dimensions by using an elliptical covariance structure. The third experiment evaluates the case where this dependency between adjacent dimensions is now dependent on which class the data is drawn from. In the fourth experiment, we return to linearly independent dimensions and consider the case where one of the distributions isn’t normally distributed, but instead uniformly distributed within a $d$-dimensional hypercube

[TABLE]

A detailed description of the distribution types and parameter setting used in each of the four experiments is presented in Table I. In addition to using the proposed method, we also estimate the Hellinger distance using two different plug-in estimators, one based on a parametric density estimator that assumes the data is normally distributed and one based on a $k$-NN density estimate of the underlying distributions. The proposed method is implemented using $\lambda=0.01$ and $k=10$, while setting $\bm{\tilde{\eta}}=[0,0.01,\dots,1]$. To calculate the $k$-NN estimate, we use the universal divergence estimation approach described in [53] and implemented in the ITE toolbox [54]. This method allows us to fix $k=10$ and still achieve an asymptotically consistent estimate of the divergence.

In each of the first three experiments, the parametric model shows the highest rate of convergence as expected, although in experiment two it is slightly outperformed at smaller sample sizes by the proposed method. In the fourth experiment, when the assumption of Gaussianity in $f_{1}$ no longer holds, the parametric solution is significantly biased and as a result, is outperformed by both non-parametric methods at higher sample sizes ($N>2000$). Relative to the $k$-NN plug-in estimator, the proposed method performs slightly worse in experiment 1, slightly better in experiment 2 and significantly better in experiments 3 and 4, with the results being relatively consistent across the various sample sizes. The improvement in performance shown in experiments 3 and 4 suggests that the proposed method offers the greatest benefit when there exists differences in the shapes of the two underlying distributions. Though the density-weighted procedure consistently outperformed the uniform method, the observed improvement was relatively minor in these experiments.

VI-B Kullback Leibler Divergence

The Kullback Leibler (KL) divergence [55], also sometimes referred to as the KL risk or relative entropy, is an asymmetric measure of divergence between two probability density functions. Using our regular notation the KL-divergence can be calculated by

[TABLE]

While the KL-divergence has the same general purpose as the Hellinger distance, that is to measure the dissimilarity between two probability density functions, it also has several key difference. Firstly since the KL-divergence is asymmetric it doesn’t technically qualify as a distance function and $d_{KL}(f_{0}||f_{1})$ isn’t necessarily equal to $d_{KL}(f_{1}||f_{0})$. This also means that $d_{KL}(f_{0}||f_{1})$ and $d_{KL}(f_{1}||f_{0})$ will have different posterior mapping functions. We can define the posterior mapping function for $d_{KL}(f_{0}||f_{1})$ as

[TABLE]

and the posterior mapping function for $d_{KL}(f_{1}||f_{0})$ as

[TABLE]

such that

[TABLE]

It is worth noting that when $p_{0}=p_{1}=0.5$

[TABLE]

thus $g_{KL}^{1}(\eta)$ is a reflection of $g_{KL}^{0}(\eta)$. One challenge presented in modeling the KL-divergence is that the posterior mapping functions are difficult to model due to discontinuities at the end points, since $g_{KL}^{0}(0)=\infty$ and $g_{KL}^{0}(1)$ is undefined though

[TABLE]

Due to their symmetry $g_{KL}^{1}$ has the same problem at the opposite endpoints. To handle this we simply select our discretized set of posteriors $\tilde{\eta}_{1},\tilde{\eta}_{2},...,\tilde{\eta}_{\tilde{N}}$ such that $0<\tilde{\eta}_{1}<\tilde{\eta}_{2}<...<\eta_{\tilde{N}}<1$. For the experiments in this Section, we set

[TABLE]

where $\epsilon=10^{-4}$. Using this modified set of discretized posteriors, we repeat the set of four experiments described in Section VI-A to evaluate the proposed methods ability to estimate the KL-divergence. The results of this experiment are displayed in Figure 8. Like the estimates of the Hellinger distance, the parametric method generally yielded the best performance in the first three experiments, but suffered from a large asymptotic bias in experiment 4. The proposed method once again significantly outperformed the $k$-NN plug-in estimator in experiments 3 and 4, however the results in experiments 1 and 2 are slightly more nuanced due to the significant difference in performance between the two optimization criteria in these trials. In both of these trials the density-weighted criteria significantly outperforms the uniform method at all sample sizes. In experiment 1 the $k$-NN plug-in estimator outperforms the regular plug-in estimator at all sample sizes, but is outperformed by the density-weighted method for $N>300$. In experiment 2 the $k$-NN plug-in estimator consistently outperforms both proposed methods, however the improvement over the density-weighted method is marginal.

VI-C $D_{p}$-Divergence

The $D_{p}$-divergence is an $f$-divergence defined by

[TABLE]

The $D_{p}$-divergence has the unique property of being directly estimable from data using minimum spanning trees [20]. Because of this property, it has been used to form non-parametric estimates of the Fisher information [40] as well as upper bounds on the Bayes error rate in a range of classification problems [20, 56]. The posterior mapping function for the $D_{p}$-divergence can be defined as

[TABLE]

which simplifies to $(2\eta-1)^{2}$ when $p_{0}=p_{1}=0.5$. We once again repeat the experiments described in Section VI-A to evaluate the proposed methods ability to estimate the $D_{p}$-divergence. This experiment provides the unique opportunity to compare the proposed method to a more traditional direct estimation procedure. The results of this experiment are displayed in Figure 9.

As in the previous experiments, the parametric estimate generally performed the best in the first three experiments, but suffered from a large asymptotic bias in experiment 4. The proposed methods perform better than the MST-based estimation in experiments 1 and 2 but worse in experiments 3 and 4. The relative performance of each method in this experiment was largely consistent across the various sample sizes, though the proposed method seems to be converging slightly faster than the MST method in experiment 4 and could possibly exceed its performance given enough samples. Unlike in the previous experiments, we found no difference in performance between the uniform optimization criteria and the density-weighted criteria in this experiment. This is due to the fact that $g_{D_{p}}(\eta)$ is a polynomial and can be perfectly represented by the proposed basis set, even when $k$ is truncated. Since we are able to achieve a solution where $\epsilon(\eta)=0$ $\forall\eta$, the density weighting has no impact on the final results.

VII Fitting Bounds on Performance

The optimization criteria in the proposed fitting routines gives us the ability to not only approximate information-theoretic functions, but to bound them as well. This is especially useful for forming bounds on the Bayes Error Rate (BER). The Bayes error rate represents the optimal classification performance that is achievable for a given pair of class distributions $f_{0}(\mathbf{x})$ and $f_{1}(\mathbf{x})$ with prior probabilities $p_{0}$ and $p_{1}$ respectively and can be calculated by

[TABLE]

In essence the BER measures the intrinsic difficulty of a particular classification problem based on the data. A thorough understanding of the BER of a particular problem can help design optimal classifiers. Because of the challenges associated with estimating the BER, much of the literature has focused on generating bounds on the BER [37, 20], which are generally formulated in terms of some measure of divergence between the two class distributions. One such bound, the well-known Bhattacharya bound, is given by [37]

[TABLE]

where

[TABLE]

While the Hellinger distance here can be estimated via any of the methods discussed in Section VI-A, parametric estimates are most common. Alternatively [20] introduced the bounds

[TABLE]

where

[TABLE]

These bounds have the advantage of being provably tighter than the Bhattacharyya bounds [20]. Furthermore since $D_{\frac{1}{2}}$ represents a particular case of the $D_{p}$-divergence, which is estimable directly from data, these bounds bypass the need for density estimation much like the approaches proposed in this paper. While these bounds are significantly tighter than the Bhattacharyya bounds, they still leave room for improvement. Avi-Itzhak proposed arbitrarily tight bounds on the BER in [39], however these bounds require density estimation to be employed in practical problems. In this section, we will use a modified version of the previously described fitting routine in order to investigate how tightly we are able to bound the BER using a linear combination of directly estimable basis functions.

Using the fitting routine described in (36) to bound the BER, requires that we define $g(\eta)$ appropriately for estimation of the BER, and constrain our fit such that

[TABLE]

We can express (53) as

[TABLE]

so g(\eta)=\min\big{[}1-\eta,\eta\big{]}. Incorporating these changes within the regularized fitting routine described in (36) yields

[TABLE]

Figure 10 compares the theoretical values of each of these bounds as a function of $\eta$. These results indicate that the proposed method offers much tighter theoretical bounds than the other two methods, however this bound is based on the asymptotic properties of the proposed basis set and doesn’t consider the limitations of a finite sample estimate.

We evaluate the finite-sample performance of this method by calculating each of the described bounds across the four experiments described in Table I. Figure 11 displays the ground truth value of the BER, along with the theoretical and estimated values for each of the three bounds (Bhattacharyya, $D_{p}$, and convex) across sample sizes ranging from 100 to 10000. The Bhattacharyya bound is calculated based on a parametric plug-in estimator which assumes both class distributions to be normally distributed. The $D_{p}$ bound is calculated from the Friedman-Rafsky test statistic using the approach described in [20]. Finally the convex bound is calculated as a linear combination of the proposed directly estimable basis functions using weights optimized according to (60). The results of this experiment are largely consistent across the four experiments, the convex method yields the tightest bound, followed by the $D_{p}$ bound, and finally the Bhattacharyya bound. Except for the Bhattacharyya bound in experiment 4, which is estimated parametrically, all of the bounds appear to converge to their asymptotic solution. While the convex bound generally offers a slightly slower convergence rate than the other two solutions, it remains tighter than the other two bounds across all sample sizes.

In order to further validate this bound we repeat one of the experiments conducted in [20] by evaluating the proposed bound along with the Mahalanobis bound, the Bhattacharyya bound, and the $D_{p}$ bound on two 8-dimensional Gaussian data sets described in [57]. The mean and standard deviations of $f_{0}$ and $f_{1}$ for the two data sets are described in Table II, and all dimensions are independent. These data sets allow us to analyze the tightness and validity of the bounds in a higher dimensional setting. For this experiment the sample size was fixed at $N=1000$ and only the empirical value of each of the bounds was evaluated. Table III displays the mean and standard deviation of each bound calculated across 500 Monte Carlo iterations for each of the two data sets. In both data sets the convex method provides the tightest bounds on the BER.

VIII Conclusion

This paper introduces a novel method for estimating density functionals which utilizes a set of directly estimable basis functions. The most appealing feature of the proposed method is its flexibility. Where previous methods of direct estimation are generally only applicable to a specific quantity, we show that the basis set can be used to generate an asymptotically consistent estimate of a broad class of density functionals, including all $f$-divergences and the Bayes error rate. We validate these findings by experimentally evaluating the proposed method’s ability to estimate three different divergences (the KL-divergence, the Hellinger distance, and the $D_{p}$-divergence) for four pairs of multivariate probability density functions. The results reveal that the proposed method performs competitively with other non-parametric divergence estimation methods, and seems to outperform them in cases where the data from the two distributions have different covariance structures or belong to different families. Additionally we demonstrate how the method can be modified to generate empirically-estimable bounds on the Bayes error rate that are much tighter than existing bounds.

Future work could focus on studying the finite-sample properties of the basis set proposed in this paper, since this represents a major source of error for the proposed methodology. An improved understanding of these properties could enable us to refine the regularization term in our optimization criteria to more accurately model each weights contribution to the estimation error or to develop ensemble methods, like those in [26], for estimating the individual basis functions. Another worthwhile future direction would be on determining an orthogonal version of the Bernstein basis expansion since this would simplify solving for the weights of the expansion.

Appendix A Proof of Theorem 1

Aspects of this proof mirror the methods used to prove the asymptotic convergence of the $k$-NN error rate in pages 62-70 of [58], please consult this text for additional details. Construction of this proof requires the introduction of an auxilary variable $\mathbf{u}$. We define the augmented data set $[\mathbf{X},\mathbf{y},\mathbf{u}]$, where $\mathbf{X}$ is as defined previously, $\mathbf{u}$ is a set of i.i.d. random variables which are uniformly distributed in $[0,1]$, and $\mathbf{y}$ is defined such that

[TABLE]

Now, let us consider an arbitrary point $\mathbf{x}^{*}$ in the support of $f_{\mathbf{x}}(\mathbf{x})$. For each instance $\mathbf{x}_{i}$, we can define an alternate label

[TABLE]

Unlike the real labels $y_{i}$, these alternate labels $y^{\prime}_{i}$ have no dependency on $\mathbf{x}_{i}$; they depend only on the fixed point $\mathbf{x}^{*}$. From these alternate labels we construct the alternate statistic

[TABLE]

Since $\Phi_{k}^{\prime}(\mathbf{x}^{*})$ is simply the sum of $k$ i.i.d. Bernoulli random variables, we can express the probability that $\Phi_{k}^{\prime}(\mathbf{x}^{*})=r$ as

[TABLE]

We can upper bound the likelihood that $\Phi^{\prime}_{k}(\mathbf{x}^{*})\neq\Phi_{k}(\mathbf{x}^{*})$ by

[TABLE]

Using the definitions of $y_{i}$ and $y_{i}^{\prime}$ above, (65) can be expressed in terms of the difference between the posterior likelihood at $\mathbf{x}^{*}$ vs. $\mathbf{x}_{i}$ as

[TABLE]

Using Lemma 5.4 in [58], we can show that

[TABLE]

as $N\rightarrow\infty$ whenever $\frac{k}{N}\rightarrow 0$. Combining (66) and (67)

[TABLE]

as $N\rightarrow\infty$ whenever $\frac{k}{N}\rightarrow 0$. Furthermore, since convergence in probability implies convergence in distribution [59],

[TABLE]

Using (69), we can simplify the expectation of $I_{r,k}(\mathbf{x})$ to

[TABLE]

where $\mathbf{x}^{\prime}$ represents a random variable independent of all instances $\mathbf{x}_{i}$ and distributed according to $f_{\mathbf{x}}(\mathbf{x}^{\prime})$. Now, let us evaluate the expression

[TABLE]

Beginning with the first term in (71)

[TABLE]

Note that the random variables $I_{r,k}(\mathbf{x}_{i})$ and $I_{r,k}(\mathbf{x}_{j})$ are conditionally independent if $B$ is true, where $B=\{\mathbf{x}_{i}\in\mathbf{X}\setminus\mathcal{N}_{k}(\mathbf{x}_{j})\cap\mathbf{x}_{j}\in\mathbf{X}\setminus\mathcal{N}_{k}(\mathbf{x}_{i})\}$. Since every instance in $\mathbf{X}\setminus\mathbf{x}_{i}$ is equally likely to be in $\mathcal{N}_{k}(\mathbf{x}_{i})$, we can bound the probability of the complement of $B$ by $P[\bar{B}]\leq 2k/N$. From this we know that as $N\rightarrow\infty$, $P[B]\rightarrow 1$ and as a result

[TABLE]

Plugging this into (72) yields

[TABLE]

Similarly the second term simplifies to

[TABLE]

and the third term clearly equals $(\rho^{*})^{2}$. Substituting these results into (71) yields

[TABLE]

Appendix B Proof of Theorem 21

Starting with the expression

[TABLE]

we first evaluate the limit with respect to $N$. Since these conditions mirror those of Theorem 1 and $g(r/k)$ is independent of $N$, we can rewrite (77) as

[TABLE]

which according to Weierstrass’ Approximation Theorem simplifies to

[TABLE]

Therefore

[TABLE]

Appendix C Proof of Theorem 3

Starting with the augmented data set defined in Appendix A, remember that $\Phi_{k}^{\prime}(\mathbf{x}^{*})$ is the sum of $k$ i.i.d. Bernoulli random variables. Therefore $\frac{1}{k}\Phi^{\prime}_{k}(\mathbf{x}^{*})$ represents the arithmetic mean of these output values, and

[TABLE]

Now, since $\Phi_{k}^{\prime}(\mathbf{x})$ must be an integer, its probability of equaling $r$ can be expressed as

[TABLE]

Taking the limit of this expression w.r.t. $k$, allows us to use (81) to form

[TABLE]

Now if we multiply each side by $k$, the right hand side takes the form of Newton’s difference quotient and can be simplified to the probability density function

[TABLE]

Since $\rho_{r,k,N}\rightarrow P[\Phi^{\prime}_{k}(\mathbf{x})=r]$ as $N$ and $k$ approach infinity in a linked manner such that $k/N\rightarrow 0$ and $\frac{r}{k}\rightarrow\eta^{*}$

[TABLE]