Ensemble estimators for multivariate entropy estimation
The problem of estimation of density functionals like entropy and mutual information has received much attention in the statistics and information theory communities. A large class of estimators of functionals of the probability density suffer from the curse of dimensionality, wherein the mean squared error (MSE) decays increasingly slowly as a function of the sample size $T$ as the dimension $d$ of the samples increases. In particular, the rate is often glacially slow of order $O(T^{-{\gamma}/{d}})$, where $\gamma>0$ is a rate parameter. Examples of such estimators include kernel density estimators, $k$-nearest neighbor ($k$-NN) density estimators, $k$-NN entropy estimators, intrinsic dimension estimators and other examples. In this paper, we propose a weighted affine combination of an ensemble of such estimators, where optimal weights can be chosen such that the weighted estimator converges at a much faster dimension invariant rate of $O(T^{-1})$. Furthermore, we show that these optimal weights can be determined by solving a convex optimization problem which can be performed offline and does not require training data. We illustrate the superior performance of our weighted estimator for two important applications: (i) estimating the Panter-Dite distortion-rate factor and (ii) estimating the Shannon entropy for testing the probability distribution of a random sample.
\section{Introduction}

Non-linear functionals of probability densities $f$ of the form $G(f) = \int g(f(x),x) f(x) dx$ arise in applications of information theory, machine learning, signal processing and statistical estimation. Important examples of such functionals include Shannon $g(f,x)=-\log(f)$ and R\'enyi $g(f,x) = f^{\alpha-1}$ entropy, and the quadratic functional $g(f,x) = f^{2}$. In these applications, the functional of interest often must be estimated empirically from sample realizations of the underlying densities. 

Functional estimation has received significant attention in the mathematical statistics community. However, estimators of functionals of multivariate probability densities $f$ suffer from mean square error (MSE) rates which typically decrease with dimension $d$ of the sample as $O(T^{-{\gamma}/{d}})$, where $T$ is the number of samples and $\gamma$ is a positive rate parameter. Examples of such estimators include kernel density estimators~\cite{kde}, $k$-nearest neighbor ($k$-NN) density estimators~\cite{fuk2}, $k$-NN entropy functional estimators~\cite{hero,kks,litt}, intrinsic dimension estimators~\cite{kks}, divergence estimators~\cite{wang}, and mutual information estimators. This slow convergence is due to the curse of dimensionality. In this paper, we introduce a simple affine combination of an ensemble of such slowly convergent estimators and show that the weights in this combination can be chosen to significantly improve the rate of MSE convergence of the weighted estimator.  In fact our ensemble averaging method can improve MSE convergence to the parametric rate $O(T^{-1})$.

Specifically, for $d$-dimensional data, it has been observed that the variance of estimators of functionals $G(f)$ decays as $O(T^{-1})$ while the bias decays as $O(T^{-1/(1+d)})$. To accelerate the slow rate of convergence of the bias in high dimensions, we propose a weighted ensemble estimator for ensembles of estimators that satisfy conditions ${\mathscr C}.1$(\ref{BE}) and ${\mathscr C}.2$(\ref{VE}) defined in Sec. II below. Optimal weights, which serve to lower the bias of the ensemble estimator to $O(T^{-1/2})$, can be determined by solving a convex optimization problem. Remarkably, this optimization problem does not involve any density-dependent parameters and can therefore be performed offline. This then ensures MSE convergence of the weighted estimator at the parametric rate of $O(T^{-1})$. 

\subsection{Related work}

When the density $f$ is $s > d/4$ times differentiable, certain estimators of functionals of the form $\int g(f(x),x) f(x) dx$, proposed by Birge and Massart~\cite{birge}, Laurent~\cite{laurent} and  Gin\'e and Mason~\cite{gini}, can achieve the parametric MSE convergence rate of $O(T^{-1})$. The key ideas in ~\cite{birge,laurent,gini} are: (i) estimation of quadratic functionals $\int f^2(x) dx$ with MSE convergence rate $O(T^{-1})$; (ii) use of kernel density estimators with kernels that satisfy the following symmetry constraints: 
\begin{equation}
\int K(x) dx =1, \hspace{0.5in} \int x^r K(x) dx = 0, 
\label{eq:symmetrickernels}
\end{equation} for $r=1,..,s$; and finally (iii) truncating the kernel density estimate so that it is bounded away from $0$. By using these ideas, the estimators proposed by ~\cite{birge,laurent,gini} are able to achieve parametric convergence rates.

In contrast, the estimators proposed in this paper require additional higher order smoothness conditions on the density, i.~e.~ the density must be $s>d$ times differentiable. However, our estimators are much simpler to implement in contrast to the estimators proposed in ~\cite{birge,laurent,gini}. In particular, the estimators in ~\cite{birge,laurent,gini} require separately estimating quadratic functionals of the form $\int f^2(x) dx$, and using truncated kernel density estimators with symmetric kernels (\ref{eq:symmetrickernels}), conditions that are not required in this paper. Our estimator is a simple affine combination of an ensemble of estimators, where the ensemble satisfies conditions ${\mathscr C}.1$ and ${\mathscr C}.2$. Such an ensemble can be trivial to implement. For instance, in this paper we show that simple uniform kernel plug-in estimators (\ref{eq:plugin}) satisfy conditions ${\mathscr C}.1$ and ${\mathscr C}.2$. 

Ensemble based methods have been previously proposed in the context of classification. For example, in both boosting~\cite{boosting} and multiple kernel learning~\cite{multipleK} algorithms, lower complexity weak learners are combined to produce classifiers with higher accuracy. Our work differs from these methods in several ways. First and foremost, our proposed method performs estimation rather than classification. An important consequence of this is that the weights we use are {\emph {data independent}}, while the weights in boosting and multiple kernel learning must be estimated from training data since they depend on the unknown distribution. 

\subsection{Organization}
The remainder of the paper is organized as follows. We formally describe the weighted ensemble estimator for a general ensemble of estimators in Section~\ref{sec:genmeth}, and specify conditions ${\mathscr C}.1$ and ${\mathscr C}.2$ on the ensemble that ensure that the ensemble estimator has a faster rate of MSE convergence. Under the assumption that conditions ${\mathscr C}.1$ and ${\mathscr C}.2$ are satisfied, we provide an MSE optimal set of weights as the solution to a  convex optimization(\ref{convexsoltheory}). Next, we shift the focus to entropy estimation in Section~\ref{sec:entropyest}, propose an ensemble of simple uniform kernel plug-in entropy estimators, and show that this ensemble satisfies conditions ${\mathscr C}.1$ and ${\mathscr C}.2$. Subsequently, we apply the ensemble estimator theory in Section~\ref{sec:genmeth} to the problem of entropy estimation using this ensemble of kernel plug-in estimators. We present simulation results in Section~\ref{sec:exp} that illustrate the superior performance of this ensemble entropy estimator in the context of (i) estimation of the Panter-Dite distortion-rate factor~\cite{riten} and (ii) testing the probability distribution of a random sample. We conclude the paper in Section~\ref{sec:con}.

\subsection*{Notation}
We will use bold face type to indicate random variables and random vectors and regular type face for constants. We denote the statistical expectation operator by the symbol ${{\mathbb{E}}}$ and the conditional expectation given random variable $\mathbf{Z}$ using the notation ${{\mathbb{E}}}_{\mathbf{Z}}$. We also define the variance operator as ${{\mathbb{V}}}[\mathbf{X}] = {{\mathbb{E}}}[(\mathbf{X}-{{\mathbb{E}}}[\mathbf{X}])^2] $ and the covariance operator as $\mathrm{Cov}[\mathbf{X},\mathbf{Y}] = {{\mathbb{E}}}[(\mathbf{X}-{{\mathbb{E}}}[\mathbf{X}])(\mathbf{Y}-{{\mathbb{E}}}[\mathbf{Y}])] $. We denote the bias of an estimator by $\mathbb{B}$.

\section{Ensemble estimators}
\label{sec:problem}
\label{sec:genmeth}

Let $\bar{l} = \{l_1,..,l_L\}$ denote a set of parameter values. For a parameterized ensemble of estimators $\{\hat{\mathbf{E}}_l\}_{l \in \bar{l}}$ of $E$, define the weighted ensemble estimator with respect to weights $w = \{w(l_1),\ldots,w(l_L)\}$ as $$\hat{\mathbf{E}}_w = \sum_{l \in \bar{l}} w(l)\hat{\mathbf{E}}_l$$ where the weights satisfy $\sum_{l \in \bar{l}} w(l) = 1$. This latter sum-to-one condition guarantees that $\hat{\mathbf{E}}_w$ is asymptotically unbiased if the component estimators $\{\hat{\mathbf{E}}_l\}_{l \in \bar{l}}$ are asymptotically unbiased. Let this ensemble of estimators $\{\hat{\mathbf{E}}_l\}_{l \in \bar{l}}$ satisfy the following two conditions:
\begin{itemize}
 \item ${\mathscr C}.1$ The bias is given by 
\begin{eqnarray}
\label{Bias_Ensemble}
\label{BE}
\mathbb{B}({\hat{\mathbf{E}}}_{l}) &=& \sum_{i \in {\cal I}} c_{i}\psi_i(l)T^{-{i/2d}} + O(1/\sqrt{T}),
\end{eqnarray}
where $c_{i}$ are constants that depend on the underlying density, ${\cal I} = \{i_1,..,i_I\}$ is a finite index set with cardinality $I<L$, $\min({\cal I}) = i_0 > 0$ and $\max({\cal I}) = i_d \leq d$, and $\psi_i(l)$ are basis functions that depend only on the estimator parameter $l$.

\item ${\mathscr C}.2$ The variance is given by 
\begin{eqnarray}
\label{Variance_Ensemble}
\label{VE}
\mathbb{V}({\hat{\mathbf{E}}}_{l}) &=& c_{v}\left({\frac{1}{T}}\right) + o\left(\frac{1}{T}\right).
\end{eqnarray}

\end{itemize}

\begin{theorem}
\label{lemma:weightedensemble}
For an ensemble of estimators $\{\hat{\mathbf{E}}_l\}_{l \in \bar{l}}$, assume that the conditions ${\mathscr C}.1$ and ${\mathscr C}.2$ hold. Then, there exists a weight vector $w_o$ such that
$${{\mathbb{E}}}[(\hat{\mathbf{E}}_{w_o} - E)^2]  = O(1/T).$$ This weight vector can be found by solving the following convex optimization problem:
\begin{equation}
\begin{aligned}
& \underset{w}{\text{minimize}}
& & ||w||_2 \\
& \text{subject to}
& & \sum_{l \in \bar{l}} w(l) = 1, \\
&&& \gamma_w(i) = \sum_{l \in \bar{l}} w(l)\psi_i(l) = 0, \; i \in {\cal I},
\end{aligned}
\label{convexsoltheory}
\end{equation}
where $\psi_i(l)$ is the basis defined in (2.1).
\end{theorem}

\begin{proof}

The bias of the ensemble estimator is given by
\begin{eqnarray}
\label{Bias_EnsembleEstimate}
\mathbb{B}({\hat{\mathbf{E}}}_{w}) &=& \sum_{i \in {\cal I}} c_{i}\gamma_w(i)T^{-{i/2d}} + O\left(\frac{||w||_1}{\sqrt{T}}\right) \nonumber \\
&=& \sum_{i \in {\cal I}} c_{i}\gamma_w(i)T^{-{i/2d}} + O\left(\frac{\sqrt{L}||w||_2}{\sqrt{T}}\right) 
\end{eqnarray}


Denote the covariance matrix of $\{\hat{\mathbf{E}}_l; l \in \bar{l}\}$ by $\Sigma_L$. Let $\bar{\Sigma}_L = \Sigma_LT$. Observe that by (\ref{Variance_Ensemble}) and the Cauchy-Schwarz inequality, the entries of $\bar{\Sigma}_L$ are $O(1)$. The variance of the weighted estimator $\hat{\mathbf{E}}_w$ can then be bounded as follows:
\begin{eqnarray}
\label{eq:var}
{{\mathbb{V}}}(\hat{\mathbf{E}}_w) &=& {{\mathbb{V}}}(\sum_{l \in \bar{l}} w_l\hat{\mathbf{E}}_l) =  w' \Sigma_L w = \frac{w' \bar{\Sigma}_L w}{T}\nonumber \\
&\leq& \frac{\lambda_{\max}(\bar{\Sigma}_L)||w||^2_2}{T} \leq \frac{trace(\bar{\Sigma}_L)||w||^2_2}{T} \leq \frac{L||w||^2_2}{T} 
\end{eqnarray}

We seek a weight vector $w$ that (i) ensures that the bias of the weighted estimator is $O(T^{-1/2})$ and (ii) has low $\ell_2$ norm $||w||_2$ in order to limit the contribution of the variance, and the higher order bias terms of the weighted estimator.  To this end, let $w_{o}$ be the solution to the convex optimization problem defined in (2.3). The solution $w_o$ is the solution of 
\begin{equation}
\begin{aligned}
& \underset{w}{\text{minimize}}
& & ||w||^2_2 \\
& \text{subject to}
& & A_0w = b, 
\end{aligned}
\nonumber
\label{convexsol2}
\end{equation}
where $A_0$ and $b$ are defined below. Let $a_0$ be the vector of ones: $[1,1. .. ,1]_{1 \times L}$; and let $a_{i}$, for each $i \in \cal{I}$ be given by $a_{i} = [\psi_i(l_1), .. ,\psi_i(l_L)]$. Define $A_0 = [a'_0, a'_{i_1}, ... , a'_{i_I}]'$, $A_1 = [a'_{i_1}, ... , a'_{i_I}]'$ and $b = [1;0;0;..;0]_{(I+1) \times 1}$. 

 Since $L > I$, the system of equations $A_0 w = b$ is guaranteed to have at least one solution (assuming linear independence of the rows $a_i$).  The minimum squared norm $\eta_L(d) := ||w_0||_2^2$ is then given by
$$\eta_L(d) = {\frac{\text{det}(A_1A'_1)}{\text{det}(A_0A'_0)}}.$$ 

Consequently, by (\ref{Bias_EnsembleEstimate}), the bias $\mathbb{B}[\hat{\mathbf{E}}_{w_o}] = O(\sqrt{L\eta_L(d)}/\sqrt{T})$. By (\ref{eq:var}), the estimator variance ${{\mathbb{V}}}[\hat{\mathbf{E}}_{w_0}] = O(L\eta_L(d)/T)$. The overall MSE is also therefore of order $O(L\eta_L(d)/T)$.

For any fixed dimension $d$ and fixed number of estimators $L>I$ in the ensemble independent of sample size $T$, the value of $\eta_L(d)$ is also independent of $T$. Stated mathematically, $L\eta_L(d) = \Theta(1)$ for any fixed dimension $d$ and fixed number of estimators $L>I$ independent of sample size $T$. This concludes the proof.

\end{proof}

In the next section, we will verify conditions $\mathscr C.1$(\ref{BE}) and $\mathscr C.2$(\ref{VE}) for plug-in estimators $\hat{\mathbf{G}}_{k}(f)$ of entropy-like functionals $G(f) = \int g(f(x),x) f(x) dx$. 


\section{Application to estimation of functionals of a density}
\label{sec:entropyest}
Our focus is the estimation of general non-linear functionals $G(f)$ of $d$-dimensional multivariate densities $f$ with known finite support ${\cal S} = [a,b]^d$, where $G(f)$ has the form
\begin{equation}
\label{eq:oracle}
G(f) = \int g(f(x),x) f(x) dx, 
\end{equation}
for some smooth function $g(f,x)$. Let ${\cal B}$ denote the boundary of ${\cal S}$. Assume that $T=N+M$ i.i.d realizations  $\{\mathbf{X}_1, \ldots, \mathbf{X}_N, \mathbf{X}_{N+1}, \ldots, \mathbf{X}_{N+M}\}$ are available from the density $f$. 

\subsection{Plug-in estimators of entropy}
\label{sec:weightedpluginest}
The truncated \emph{uniform kernel} density estimator is defined below. For any positive real number $k \leq M$, define the distance $d_k$ to be: $d_k = (k/M)^{1/d}$. Define the truncated kernel region for each $X \in {\cal S}$ to be $S_k(X) = \{Y \in {\cal S} : ||X-Y||_\infty \leq d_k/2\}$, and the volume of the truncated uniform kernel to be $V_k(X) = \int_{S_k(X)} dz$. Note that when the smallest distance from $X$ to ${\cal B}$ is greater than $d_k/2$, $V_k(X) = d_k^d = k/M$. Let $\mathbf{l}_k(X)$ denote the number of samples falling in $S_k(X)$: $\mathbf{l}_k(X) = \sum_{i=1}^{M} 1_{\{\mathbf{X}_i \in S_k(X)\}}$. The truncated {uniform kernel} density estimator is defined as
\begin{equation}
  \hat{\mathbf{f}}_{k}(X) = \frac{\mathbf{l}_k(X)}{MV_k(X)}.
\end{equation}

The plug-in estimator of the density functional is constructed using a data splitting approach as follows. The data is randomly subdivided into two parts $\{\mathbf{X}_1, \ldots, \mathbf{X}_N\}$ and $\{\mathbf{X}_{N+1}, \ldots, \mathbf{X}_{N+M}\}$ of $N$ and $M$ points respectively. In the first stage, we form the kernel density estimate ${\hat{\mathbf{f}}_k}$ at the $N$ points $\{\mathbf{X}_1, \ldots, \mathbf{X}_N\}$ using the $M$ realizations $\{\mathbf{X}_{N+1}, \ldots, \mathbf{X}_{N+M}\}$. Subsequently, we use the $N$ samples $\{\mathbf{X}_1, \ldots, \mathbf{X}_N\}$ to approximate the functional $G(f)$ and obtain the plug-in estimator:
\begin{eqnarray}
\label{eq:plugin}
  \hat{\mathbf{G}}_k &=& \frac{1}{N}\sum_{i=1}^N g({\hat{\mathbf{f}}{_k}(\mathbf{X}_i)},\mathbf{X}_i). 
\end{eqnarray}
Also define a standard kernel density estimator $\tilde{\mathbf{f}}_k$, which is identical to $\hat{\mathbf{f}}_k$ except that the volume $V_k(X)$ is always set to the untruncated value $V_k(X) = k/M$. Define
\begin{eqnarray}
\label{eq:plugin2}
  \tilde{\mathbf{G}}_k &=& \frac{1}{N}\sum_{i=1}^N g({\tilde{\mathbf{f}}{_k}(\mathbf{X}_i)},\mathbf{X}_i). 
\end{eqnarray}
The estimator $\tilde{\mathbf{G}}_k$ is identical to the estimator of Gy{\"{o}}rfi and van der Meulen~\cite{gyrofi}. Observe that the implementation of $\tilde{\mathbf{G}}_k$, unlike $\hat{\mathbf{G}}_k$, does not require knowledge about the support of the density.

 

\subsubsection{Assumptions}
\label{sec:assump}  
We make a number of technical assumptions that will allow us to obtain tight MSE convergence rates for the kernel density estimators defined above. $({\cal {A}}.0)$ : Assume that $k = k_0M^\beta$ for some rate constant $0<\beta<1$, and assume that $M$, $N$ and $T$ are linearly related through the proportionality constant $\alpha_{frac}$ with: $0 < \alpha_{frac} < 1$, $M = \alpha_{frac}T$ and $N = (1-\alpha_{frac})T$. $({\cal {A}}.1)$ : Let the density $f$ be uniformly bounded away from $0$ and upper bounded on the set ${\cal S}$, i.e., there exist constants $\epsilon_0$, $\epsilon_\infty$ such that $0 < \epsilon_0 \leq f(x) \leq \epsilon_\infty < \infty$ $\forall x \in {\cal S}$. $({\cal {A}}.2)$: Assume that the density $f$ has continuous partial derivatives of order $d$ in the interior of the set ${\cal S}$, and that these derivatives are upper bounded. $({\cal {A}}.3)$: Assume that the function $g(f,x)$ has $\max\{\lambda,d\}$ partial derivatives w.r.t. the argument $f$, where $\lambda$ satisfies the condition  $\lambda \beta>1$. Denote the $n$-th partial derivative of $g(f,x)$ wrt $x$ by $g^{(n)}(f,x)$.  $({\cal {A}}.4)$: Assume that the absolute value of the functional $g(f,x)$ and its partial derivatives are strictly upper bounded in the range $\epsilon_0 \leq f \leq \epsilon_\infty$ for all $x$. $({\cal {A}}.5)$: Let $\epsilon \in (0,1)$ and $\delta \in (2/3,1)$. Let ${\cal C}(M)$ be a positive function satisfying the condition $ {\cal C}(M) = \Theta(\exp(-M^{\beta(1-\delta)}))$. For some fixed $0 < \epsilon< 1$, define $p_l = (1-\epsilon)\epsilon_0$ and $p_u = (1+\epsilon)\epsilon_\infty $. Assume that the conditions
$$ (i) \sup_{x}|h(0,x)|  < G_1 < \infty, $$
$$ (ii) \sup_{f \in (p_l,p_u),x}|h(f,x)|  < G_2 < \infty, $$
$$(iii) \sup_{f \in (1/k,p_u),x }|h(f,x)|{\cal C}(M)  < G_3 < \infty \hspace{0.15in} \forall M,$$ 
$$(iv)\sup_{f \in (p_l,2^dM/k),x} |h(f,x)|{\cal C}(M)  < G_4 < \infty \hspace{0.15in} \forall M,$$ are satisfied by 
$h(f,x) = g(f,x), g^{(3)}(f,x)$ and $g^{(\lambda)}(f,x)$, for some constants $G_1$, $G_2$, $G_3$ and $G_4$.

These assumptions are comparable to other rigorous treatments of entropy estimation. The assumption $({\cal {A}}.0)$ is equivalent to choosing the bandwidth of the kernel to be a fractional power of the sample size~\cite{raykar}. The rest of the above assumptions can be divided into two categories: (i) assumptions on the density $f$, and (ii) assumptions on the functional $g$. The assumptions on the smoothness, boundedness away from $0$ and $\infty$ of the density $f$ are similar to the assumptions made by other estimators of entropy as listed in Section II,~\cite{beir}. The assumptions on the functional $g$ ensure that $g$ is sufficiently smooth and that the estimator is bounded. These assumptions on the functional are readily satisfied by the common functionals that are of interest in literature: Shannon $g(f,x) = - \log(f)I(f>0) + I(f=0)$ and R\'enyi $g(f,x) = f^{\alpha-1}I(f>0) + I(f=0)$ entropy, where $I(.)$ is the indicator function, and the quadratic functional $g(f,x) = f^{2}$.


\subsubsection{Analysis of MSE}
\label{sec:mseanal}
Under the assumptions stated above, we have shown the following in the Appendix:
\begin{theorem}
\label{knnbiasH}
The biases of the plug-in estimators $\hat{\mathbf{G}}_{k}, \tilde{\mathbf{G}}_{k}$ are given by
\begin{eqnarray}
\label{Biasu}
\mathbb{B}(\hat{\mathbf{G}}_{k}) &=& \sum_{i=1}^{d} c_{1,i}\left({\frac{k}{M}}\right)^{i/d} + \frac{c_{2}}{k} \nonumber + o\left(\frac{1}{k} + \frac{k}{M}\right) \nonumber \\
\mathbb{B}(\tilde{\mathbf{G}}_{k}) &=& c_{1}\left({\frac{k}{M}}\right)^{1/d} + \frac{c_{2}}{k} \nonumber + o\left(\frac{1}{k} + \frac{k}{M}\right), \nonumber
\end{eqnarray}
where $c_{1,i}$, $c_1$ and $c_2$ are constants that depend on $g$ and $f$
\end{theorem}

\begin{theorem}
\label{knnvarH}
The variances of the plug-in estimators $\hat{\mathbf{G}}_{k}, \tilde{\mathbf{G}}_{k}$ are identical up to leading terms, and are given by
\begin{eqnarray}
\label{Variance}
{{\mathbb{V}}}(\hat{\mathbf{G}}_{k}) &=& c_4\left(\frac{1}{N}\right)+ c_5\left(\frac{1}{M}\right) + o\left(\frac{1}{M} + \frac{1}{N}\right) \nonumber \\
{{\mathbb{V}}}(\tilde{\mathbf{G}}_{k}) &=& c_4\left(\frac{1}{N}\right)+ c_5\left(\frac{1}{M}\right) + o\left(\frac{1}{M} + \frac{1}{N}\right), \nonumber
\end{eqnarray}
where $c_4$ and $c_5$ are constants that depend on $g$ and $f$.
\end{theorem}


\subsubsection{Optimal MSE rate}

From Theorem \ref{knnbiasH}, observe that the conditions $k \to \infty$ and $k/M \to 0$ are necessary for the estimators $\hat{\mathbf{G}}_{k}$ and $\tilde{\mathbf{G}}_{k}$ to be unbiased. Likewise from Theorem \ref{knnvarH}, the conditions $N \to \infty$ and $M \to \infty$ are necessary for the variance of the estimator to converge to $0$. Below, we optimize the choice of bandwidth $k$ for minimum MSE, and also show that the optimal MSE rate is invariant to the choice of $\alpha_{frac}$.

\paragraph{Optimal choice of $k$}
Minimizing the MSE over $k$ is equivalent to minimizing the square of the bias over $k$. The optimal choice of $k$ is given by
\begin{eqnarray}
\label{kopt}
k_{opt} &=& \Theta({M^{{1}/{1+d}}}),
\end{eqnarray}
and the bias evaluated at $k_{opt}$ is  $\Theta({M^{{-1}/{1+d}}})$.

\paragraph{Choice of $\alpha_{frac}$}
Observe that the MSE of $\hat{\mathbf{G}}_{k}$ and $\tilde{\mathbf{G}}_{k}$ are dominated by the squared bias $(\Theta(M^{-2/(1+d)}))$ as contrasted to the variance $(\Theta(1/N+1/M))$. This implies that the asymptotic MSE rate of convergence is invariant to the selected proportionality constant  $\alpha_{frac}$.

In view of (a) and (b) above, the optimal MSE for the estimators $\hat{\mathbf{G}}{_k}$ and $\tilde{\mathbf{G}}{_k}$ is therefore achieved for the choice of $k=\Theta(M^{1/(1+d)})$, and is given by $\Theta(T^{-2/(1+d)})$. Our goal is to reduce the estimator MSE to $O(T^{-1})$. We do so by applying the method of weighted ensembles described in Section~\ref{sec:genmeth}.

\subsection{Weighted ensemble entropy estimator}
For a positive integer $L > I = d-1$, choose $\bar{l} = \{l_1, \ldots, l_L\}$ to be positive real numbers. Define the mapping $k(l) = l\sqrt{M}$ and let $\bar{k} =  \{k(l); l \in \bar{l}\}$. Define the weighted ensemble estimator 
\begin{equation}
\hat{\mathbf{G}}_w  = \sum_{l \in \bar{l}} w(l)\hat{\mathbf{G}}_{k(l)}.
\label{eq:weightedensemble}
\end{equation}
From Theorems \ref{knnbiasH} and \ref{knnvarH}, we see that the biases of the ensemble of estimators $\{\hat{\mathbf{G}}_{k(l)}; l \in \bar{l}\}$ satisfy ${\mathscr C}.1$(\ref{BE}) when we set  $\psi_i(l) = l^{i/d}$ and ${\cal I} = \{1,..,d-1\}$. Furthermore, the general form of the variance of $\hat{\mathbf{G}}_{k(l)}$ follows ${\mathscr C}.2$(\ref{VE}) because $N, M = \Theta(T)$. This implies that we can use the weighted ensemble estimator $\hat{\mathbf{G}}_w$ to estimate entropy at $O(L\eta_L(d)/T)$ convergence rate by setting $w$ equal to the optimal weight $w_o$ given by (\ref{convexsoltheory}).



\section{Experiments}
\label{sec:exp}

We illustrate the superior performance of the proposed weighted ensemble estimator for two applications: (i) estimation of the Panter-Dite rate distortion factor, and (ii) estimation of entropy to test for randomness of a random sample. 

For finite $T$ direct use of Theorem 1 can lead to excessively high variance. This is because forcing the condition (2.3) that $\gamma_w(i)=0$ is too strong and, in fact, not necessary. The careful reader may notice that to obtain $O(T^{-1}) $ MSE convergence  rate in Theorem 1 it is sufficient that $\gamma_w(i)$ be of order $O(T^{-1/2+i/2d})$. Therefore, in practice we determine the optimal weights according to the optimization:

\begin{equation}
\begin{split}
\min_w \quad &\epsilon\\
\text{subject to} \quad &\gamma_w(0) = 1,\\
&\lvert \gamma_w(i) T^{1/2 - i/2d} \rvert \leq \epsilon, \quad i \in \mathcal{I},\\
&\lVert w \rVert_2^2 \leq \eta.
\end{split}
\label{convexsol}
\end{equation}

The optimization \eqref{convexsol} is also  convex. Note that, as contrasted to \eqref{convexsoltheory}, the norm of the weight vector $w$ is bounded instead of being minimized.  By relaxing the constraints $\gamma_w(i)=0$ in \eqref{convexsoltheory} to the softer constraints in \eqref{convexsol}, the upper bound $\eta$ on $\lVert w \rVert_2^2$ can be reduced from the value $\eta_L(d)$ obtained by solving \eqref{convexsoltheory}.  This results in a more favorable trade-off between bias and variance for moderate sample sizes.  In our experiments, we find that setting $\eta = 3d$ yields good MSE performance.  Note that as $T \to \infty$, we must have $\gamma_w(i) \to 0$ for $i \in \mathcal{I}$  in order to keep $\epsilon$ finite, thus recovering the strict constraints in \eqref{convexsoltheory}.



For fixed sample size $T$ and dimension $d$, observe that increasing $L$ increases the number of degrees of freedom in the convex problem \eqref{convexsol}, and therefore will result in a smaller value of $\epsilon$ and in turn improved estimator performance.   In our simulations, we choose $\bar{l}$ to be $L=50$ equally spaced values between $0.3$ and $3$, ie the $l_i$ are uniformly spaced as $$l_i = \frac{x}{a} + \frac{(a-1)ix}{aL}; i=1,..,L,$$ with scale and range parameters $a=10$ and $x=3$ respectively. We limit $L$ to 50 because we find that the gains beyond $L=50$ are negligible. The reason for this diminishing return is a direct result of the increasing similarity among the entries in $\bar{l}$, which translates to increasingly similar basis functions $\psi_i(l) = l^{i/d}$. 



\subsection{Panter-Dite factor estimation}

\begin{figure}[!t]
\centering
\subfigure[\small{Variation of MSE of Panter-Dite factor estimates as a function of sample size $T$. From the figure, we see that the proposed weighted estimator has the fastest MSE rate of convergence wrt sample size $T$ ($d=6$).}]{
\includegraphics[scale=.25]{MSEvsT.png}
\label{a-compare}
}
\subfigure[\small{Variation of MSE of Panter-Dite factor estimates as a function of dimension $d$. From the figure, we see that the MSE of the proposed weighted estimator has the slowest rate of growth with increasing dimension $d$ ($T=3000$).}]{
\includegraphics[scale=.25]{MSEvsd.png}
\label{ad-compare}
}
\caption{Variation of MSE of Panter-Dite factor estimates using standard kernel plug-in estimator~[14], truncated kernel plug-in estimator~(3.3), histogram plug-in estimator[17], $k$-NN estimator~[20], entropic graph estimator~[18] and the weighted ensemble estimator~(3.6). }
\end{figure}

For a $d$-dimensional source with underlying density $f$, the Panter-Dite distortion-rate function~\cite{riten} for a $q$-dimensional vector quantizer with $n$ levels of quantization is given by $ \delta(n)= n^{-2/q} \int f^{q/(q+2)}(x) dx. $ The Panter-Dite factor corresponds to the functional $G(f)$ with $g(f,x) = n^{-2/q} f^{-2/(q+2)}I(f>0) + I(f=0)$. The Panter-Dite factor is directly related to the R\'{e}nyi $\alpha$-entropy, for which several other estimators have been proposed~\cite{gyrofi0,heroJ,pal,leo2}. 

In our simulations we compare six different choices of functional estimators - the three estimators previously introduced: (i) the standard kernel plug-in estimator $\tilde{\mathbf{G}}_k$, (ii) the boundary truncated plug-in estimator $\hat{\mathbf{G}}_k$ and (iii) the weighted estimator $\hat{\mathbf{G}}_w$ with optimal weight $w = w^*$ given by (\ref{convexsol}), and in addition the following popular entropy estimators: (iv) histogram plug-in estimator~\cite{gyrofi0}, (v) $k$-nearest neighbor ($k$-NN) entropy estimator~\cite{leo2} and (vi) entropic $k$-NN graph estimator~\cite{heroJ,pal}. For both $\tilde{\mathbf{G}}_k$ and $\hat{\mathbf{G}}_k$, we select the bandwidth parameter $k$ as a function of $M$ according to the optimal proportionality $k=M^{1/(1+d)}$ and $N=M=T/2$. 

We choose $f$ to be the $d$ dimensional mixture density $f(a,b,p,d) = pf_\beta(a,b,d) + (1-p)f_u(d)$; where $d=6$, $f_\beta(a,b,d)$ is a $d$-dimensional Beta density with parameters $a=6,b=6$, $f_u(d)$ is a $d$-dimensional uniform density and the mixing ratio $p$ is $0.8$. The reason we choose the beta-uniform mixture for our experiments is because it trivially satisfies all the assumptions on the density $f$ listed in Section 3.1, including the assumptions of finite support and strict boundedness away from 0 on the support. The true value of the Panter-Dite factor $\delta(n)$ for the beta-uniform mixture is calculated using numerical integration methods via the 'Mathematica' software (http://www.wolfram.com/mathematica/). Numerical integration is used because evaluating the entropy in closed form for the beta-uniform mixture is not tractable. 

The MSE values for each of the six estimators are calculated by averaging the squared error $[\hat{\delta}_i(n) - \delta(n)]^2$, $i=1,..,m$ over $m = 1000$ Monte-Carlo trials, where each $\hat{\delta}_i(n)$ corresponds to an independent instance of the estimator.





\subsubsection{Variation of MSE with sample size $T$}
The MSE results of the different estimators are shown in Fig.~\ref{a-compare} as a function of sample size $T$, for fixed dimension $d=6$. It is clear from the figure that the proposed ensemble estimator $\hat{\mathbf{G}}_w$ has significantly faster rate of convergence while the MSE of the rest of the estimators, including the truncated kernel plug-in estimator, have similar, slow rates of convergence. It is therefore clear that the proposed optimal ensemble averaging significantly accelerates the MSE convergence rate.



\subsubsection{Variation of MSE with dimension $d$}

For fixed sample size $T$ and fixed number of estimators $L$, it can be seen that $\epsilon$ increases monotonically with $d$. This follows from the fact that the number of constraints in the convex problem \ref{convexsol} is equal to $d+1$ and each of the basis functions $\psi_i(l) = l^{i/d}$ monotonically approaches $1$ as $d$ grows, . This in turn implies that for a fixed sample size $T$ and number of estimators $L$, the overall MSE of the ensemble estimator should increase monotonically with the dimension $d$.

The MSE results of the different estimators are shown in Fig.~\ref{ad-compare} as a function of dimension $d$, for fixed sample size $T=3000$. For the standard kernel plug-in estimator and truncated kernel plug-in estimator, the MSE increases  rapidly with $d$ as expected. The MSE of the histogram and $k$-NN estimators increase at a similar rate, indicating that these estimators suffer from the curse of dimensionality as well. On the other hand, the MSE of the weighted estimator also increases with the dimension as predicted, but at a slower rate. Also observe that the MSE of the weighted estimator is smaller than the MSE of the other estimators for all dimensions $d>3$.

\subsection{Distribution testing}


\begin{figure}[!t]
\centering
\subfigure[\small{Entropy estimates for random samples corresponding to hypothesis $H_0$ (experiments 1-500) and $H_1$ (experiments 501-1000).}]{
\includegraphics[scale=.30]{responsemod.png}
\label{c-compare}
}
\subfigure[\small{Histogram envelopes of entropy estimates for random samples corresponding to hypothesis $H_0$ (blue) and $H_1$ (red).}]{
\includegraphics[scale=.30]{histmod.png}
\label{d-compare}
}
\caption{Entropy estimates using standard kernel plug-in estimator, truncated kernel plug-in estimator and the weighted estimator, for random samples corresponding to hypothesis $H_0$ and $H_1$. The weighted estimator provides better discrimination ability by suppressing the bias, at the cost of some additional variance.}
\end{figure}


In this section, we illustrate the weighted ensemble estimator for non-parametric estimation of Shannon differential entropy.  The Shannon differential entropy is given by $G(f)$ where $g(f,x) = - \log(f)I(f>0) + I(f=0)$. The improved accuracy of the weighted ensemble estimator is demonstrated in the context of hypothesis testing using estimated entropy as a statistic to test for the underlying probability distribution of a random sample. Specifically, the samples under the null and alternate hypotheses $H_0$ and $H_1$ are drawn from the probability distribution $f(a,b,p,d)$, described in Section IV.A, with fixed $d=6$, $p=0.75$ and two sets of values of $a,b$ under the null and alternate hypothesis, $H_0: a=a_0,b=b_0$ versus $H_1: a=a_1,b=b_1$.

First, we fix $a_0=b_0=6$ and $a_1=b_1=5$. The density under the null hypothesis $f(6,6,0.75,6)$ has greater curvature relative to $f(5,5,0.75,6)$ and therefore has smaller entropy.   Five hundred (500) experiments are performed under each hypothesis with each experiment consisting of 1000 samples drawn from the corresponding distribution. The true entropy and  estimates $\tilde{\mathbf{G}}_k$, $\hat{\mathbf{G}}_k$ and $\hat{\mathbf{G}}_w$ obtained from each instance of $10^3$ samples are shown in Fig.~\ref{c-compare} for the 1000 experiments. This figure suggests that the  ensemble weighted estimator provides better discrimination ability by suppressing the bias, at the cost of some additional variance.

To demonstrate that the weighted estimator provides better discrimination, we plot the histogram envelope of the entropy estimates using standard kernel plug-in estimator, truncated kernel plug-in estimator and the weighted estimator for the cases corresponding to the hypothesis $H_0$ (color coded {blue}) and $H_1$ (color coded {red}) in Fig.~\ref{d-compare}. Furthermore, we quantitatively measure the discriminative ability of the different estimators using the deflection statistic $ds = {|\mu_1-\mu_0|}/{\sqrt{\sigma_0^2+\sigma_1^2}},$ where $\mu_0$ and $\sigma_0$ (respectively  $\mu_1$ and $\sigma_1$) are the sample mean and standard deviation of the entropy estimates. The deflection statistic was found to be $1.49$, $1.60$ and $1.89$ for the standard kernel plug-in estimator, truncated kernel plug-in estimator and the weighted estimator respectively. The receiver operating curves (ROC) for this entropy-based test using the three different estimators are shown in Fig.~\ref{b-compare}. The corresponding areas under the ROC curves (AUC) are given by $0.9271$, $0.9459$ and $0.9619$. 

\begin{figure}[!t]
\centering
\subfigure[\small{ROC curves corresponding to entropy estimates obtained using standard and truncated kernel plug-in estimators and the weighted estimator. The corresponding AUC are given by $0.9271$, $0.9459$ and $0.9619$.}]{
\includegraphics[scale=.30]{ROC1000.pdf}
\label{b-compare}
}
\subfigure[\small{Variation of AUC curves vs $\delta (= a_0-a_1, b_0-b_1)$ corresponding to Neyman-Pearson omniscient test, entropy estimates using the standard and truncated kernel plug-in estimators and the weighted estimator.}]{
\includegraphics[scale=.30]{auc10000.pdf}
\label{e-compare}
}
\caption{Comparison of performance in terms of ROC for the distribution testing problem. The weighted estimator uniformly outperforms the individual plug-in estimators.}
\end{figure}

In our final experiment, we fix $a_0=b_0=10$ and set $a_1=b_1=10-\delta$,   perform 500 experiments each under the null and alternate hypotheses with samples of size 5000, and plot the AUC as $\delta$ varies from $0$ to $1$ in Fig.~\ref{e-compare}. For comparison, we also plot the AUC for the Neyman-Pearson likelihood ratio test. The Neyman-Pearson likelihood ratio test, unlike the Shannon entropy based tests, is an omniscient test that assumes knowledge of both the underlying beta-uniform mixture parametric model of the density and the parameter values $a_0$, $b_0$ and $a_1$, $b_1$ under the null and alternate hypothesis respectively. Figure 4 shows that the weighted estimator {\emph {uniformly and significantly}} outperforms the individual plug-in estimators and comes closest to the performance of the omniscient Neyman-Pearson likelihood test. The relatively superior performance of the Neyman-Pearson likelihood test is due to the fact that the weighted estimator is a nonparametric estimator that has marginally higher variance (proportional to $||w^*||_2^2$) as compared to the underlying parametric model for which the Neyman-Pearson test statistic provides the most powerful test. 

\section{Conclusions}
\label{sec:con}
We have proposed a new estimator of functionals of a multivariate density based on weighted ensembles of kernel density estimators.  For ensembles of estimators that satisfy general conditions on bias and variance as specified by ${\mathscr C}.1$(\ref{BE}) and ${\mathscr C}.2$(\ref{VE}) respectively,  the weight optimized ensemble estimator has parametric $O(T^{-1})$ MSE convergence rate that can be much faster than the rate of convergence of any of the individual estimators in the ensemble. The optimal weights are determined as a solution to a convex optimization problem that can be performed offline and does not require training data. We illustrated this estimator for  uniform kernel plug-in estimators and demonstrated the superior performance of the weighted ensemble entropy estimator for (i)  estimation of the Panter-Dite factor and (ii) non-parametric hypothesis testing. 

Several extensions of the framework of this paper are being pursued: (i) using $k$-nearest neighbor ($k$-NN) estimators in place of kernel estimators; (ii) extending the framework to the case where support ${\cal S}$ {\emph {is not known}}, but for which conditions ${\mathscr C}.1$ and ${\mathscr C}.2$ hold; (iii) using ensemble estimators for estimation of other functionals of probability densities including divergence, mutual information and intrinsic dimension; and (iv) using an $l_1$ norm $\|w\|_1$in place of the $l_2$ norm $\|w\|_2$ in  the weight optimization algorithm (2.3) so as to introduce sparsity into the weighted ensemble. 

\section*{Acknowledgement}
This work was partially supported by (i) ARO grant W911NF-12-1-0443 and (ii) NIH grant 2P01CA087634-06A2.

\begin{figure}[h]
\centering
\includegraphics[width=6.2in]{sss.pdf}
\caption{Illustration for the proof of Lemma~\ref{biaslemma}.}
\label{i-compare}
\end{figure}


2013-03-05T02:01:58.000
1203
1203.5829
en
https://arxiv.org/abs/1203.5829
The problem of estimation of density functionals like entropy and mutual information has received much attention in the statistics and information theory communities. A large class of estimators of functionals of the probability density suffer from the curse of dimensionality, wherein the mean squared error (MSE) decays increasingly slowly as a function of the sample size $T$ as the dimension $d$ of the samples increases. In particular, the rate is often glacially slow of order $O(T^{-{\gamma}/{d}})$, where $\gamma>0$ is a rate parameter. Examples of such estimators include kernel density estimators, $k$-nearest neighbor ($k$-NN) density estimators, $k$-NN entropy estimators, intrinsic dimension estimators and other examples. In this paper, we propose a weighted affine combination of an ensemble of such estimators, where optimal weights can be chosen such that the weighted estimator converges at a much faster dimension invariant rate of $O(T^{-1})$. Furthermore, we show that these optimal weights can be determined by solving a convex optimization problem which can be performed offline and does not require training data. We illustrate the superior performance of our weighted estimator for two important applications: (i) estimating the Panter-Dite distortion-rate factor and (ii) estimating the Shannon entropy for testing the probability distribution of a random sample.
Statistics Theory (math.ST); Methodology (stat.ME)
Ensemble estimators for multivariate entropy estimation
Qwen/Qwen-72B
1. YES
2. YES