Empirical Squared Hellinger Distance Estimator and Generalizations to a Family of α-Divergence Estimators

Abstract

We present an empirical estimator for the squared Hellinger distance between two continuous distributions, which almost surely converges. We show that the divergence estimation problem can be solved directly using the empirical CDF and does not need the intermediate step of estimating the densities. We illustrate the proposed estimator on several one-dimensional probability distributions. Finally, we extend the estimator to a family of estimators for the family of $\alpha$-divergences, which almost surely converge as well, and discuss the uniqueness of this result. We demonstrate applications of the proposed Hellinger affinity estimators to approximately bounding the Neyman–Pearson regions.

1. Introduction

We present an empirical estimator for the squared Hellinger distance between two continuous distributions. The work is a direct extension of Perez-Cruz [1] where they provided an empirical KL divergence estimator. Their work is built upon previous works on divergence estimators such as [2,3,4,5,6]. Similar to their estimator, given two samples from two distributions, our estimator does not need to estimate the probability density functions explicitly before estimating the squared Hellinger distance between the two distributions, which makes it simple and fast. We show that the estimator converges to the true squared Hellinger distance almost surely as the sample size increases. We then extend our estimator to the family of $\alpha$-divergences, to which the squared Hellinger distance belongs. For each of the estimators, we can obtain a reverse estimator using the other direction of the two data samples, and we can also obtain a symmetric estimator by averaging the two one-sided estimators. We present several numerical examples to show the convergence of our estimators. Our newly proposed estimators can be used efficiently to approximate the adjacency of two data samples, leading to various applications in many fields of research.

2. Preliminaries on Divergences between Probability Distributions

Recall that the definition of squared Hellinger distance [7] is (for univariate continuous distributions):

$H^2(P,Q)=\frac{1}{2}{\int }_x{\left(\sqrt{p\left(x\right)}-\sqrt{q\left(x\right)}\right)}^{2}dx.$

It is symmetric and always bounded between 0 and 1.

Additionally, recall the definition of Kullback–Leibler divergence [8] is (for univariate continuous distributions):

$D_{KL}\left(P\right|\left|Q\right)={\int }_{x}p\left(x\right)log\frac{p\left(x\right)}{q\left(x\right)}dx.$

KL divergence and the squared Hellinger distance both belong to a family of f-divergences, which are central to information theory and statistics. Compared with KL divergence, the squared Hellinger distance is symmetric, and Hellinger distance forms a bounded metric between 0 and 1 on the space of probability distributions. Hellinger distance is related to total variation distance as:

$H^2(P,Q)\le TVD(P,Q)\le \sqrt{2}H(P,Q),$

where total variation distance (TVD) is defined as:

$TVD(P,Q)=\frac{1}{2}{\int }_x|p\left(x\right)-q\left(x\right)|dx.$

The squared Hellinger distance is also closely related to KL divergence and can be bounded by:

$2H^2(P,Q)\le D_{KL}\left(P\right|\left|Q\right).$

It is also a known result that KL divergence is stronger than Hellinger distance in the sense that convergence in KL divergence implies convergence in Hellinger distance, which further implies convergence in total variation distances. Therefore, Hellinger distance represents a middle ground between KL divergence and total variation distance; it is weaker than KL divergence but stronger than total variation distance in terms of convergence. As shown before, Hellinger distance has close connections to the total variation distance, which is exactly what inference depends on (KL divergence does not admit a useful lower bound on the TVD). It has another attractive property compared with KL divergence, which is the fact that the squared Hellinger distance is always bounded between zero and one for probability distributions that may or may not have the same support, whereas the KL divergence becomes infinite for probability distributions of different supports. In fact, KL divergence can be unbounded for probability distributions supported on the real line. For example, consider P to be the standard Cauchy distribution and Q to be the standard normal distribution, then $D_{KL}\left(P\right|\left|Q\right)$ diverges to infinity. Hence, an empirical estimator for KL divergence does not provide meaningful estimates in such a case, while the squared Hellinger distance is always bounded. Due to these desirable properties, we focus mainly on the squared Hellinger distance in this work. The squared Hellinger distance is a member of the family of $\alpha$-divergences (up to a scaling factor), which are defined in Cichocki and Amari [9] for $\alpha \in (0,1)$ as,

$D_A^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha }+\frac{1}{1-\alpha }-\frac{1}{\alpha (1-\alpha )}{\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx.$

The $\alpha$-divergence can also be related to TVD through the following inequalities, similar to squared Hellinger distance up to a scaling factor (see for example [10,11,12]),

$\alpha (1-\alpha )D_A^{\alpha }\left(P\right|\left|Q\right)\le TVD(P,Q)\le \sqrt{\frac{D_A^{\alpha }\left(P\right|\left|Q\right)}{2}}.$

3. Review of Empirical Sample-Based Kullback–Leibler Divergence Estimator of Continuous Distributions

Let $\mathcal{X}={\left\{x_i\right\}}_{i=1}^n,{\mathcal{X}}^{\prime }={\left\{x_j^{\prime }\right\}}_{j=1}^m$ be iid samples from P and Q in increasing order. Recall that the definition of the empirical CDFs of P and Q are, respectively,

$P_e\left(x\right)=\frac{1}{n}\sum _{i=1}^{n}U(x-x_i);Q_e\left(x\right)=\frac{1}{m}\sum _{j=1}^{m}U(x-x_j^{\prime }),$

where $U\left(x\right)$ is a unit-step function with $U\left(0\right)=0.5$. The continuous piece-wise linear interpolation of the empirical CDF of P is denoted as $P_c\left(x\right)$. It is zero for any point smaller than a joint lower bound $x_0<inf\{\mathcal{X},{\mathcal{X}}^{\prime }\}$ of the data samples from $P,Q$, and is one for anything greater than or equal to a joint upper bound $x_{n+1}>sup\{\mathcal{X},{\mathcal{X}}^{\prime }\}$ of the data samples from $P,Q$; everywhere in the middle, it is defined as:

$P_c\left(x\right)=a_{i}x+b_i,x_{i-1}<x<x_i,$

where coefficients $a_i,b_i$ are set so that $P_c\left(x\right)$ matches the values of $P_e\left(x\right)$ at the sampled values $x_i,i=1,\dots ,n$. Similarly, we can define the interpolated empirical CDF for Q, denoted as $Q_c\left(x\right)$. These empirical CDFs converge uniformly and are independent of the distribution of their CDFs.

Perez-Cruz [1] proposed an empirical KL estimator:

$\widehat{D}\left(P\right|\left|Q\right)=\frac{1}{n}\sum _{i=1}^{n}log\frac{\delta P_c\left(x_i\right)}{\delta Q_c\left(x_i\right)},$

where $\delta P_c\left(x_i\right)=(P_c\left(x_i\right)-P_c(x_i-ϵ))/ϵ$ for any $ϵ<{min}_i\{x_i-x_{i-1}\}$ denotes the left slope of $P_c$ at $x_i$ and $\delta Q_c\left(x_i\right)$ denotes the left slope of $Q_c$ at $x_i$. Here, $n=\left|\mathcal{X}\right|$ and $x_i$ are the samples from the P distribution. Ref. [1] showed that $\widehat{D}\left(P\right|\left|Q\right)-1\to D\left(P\right|\left|Q\right)$, almost surely. For this 1-D data setting, an experiment showing the convergence of their estimator is shown in Figure 1 where we plotted estimated values against increasing sample sizes, where $P,Q$ are taken to be normal distributions $N(0,1)$ and $N(1,1)$ respectively.

It is worth mentioning that the major innovation and strength of these types of empirical estimators is the fact that there are no convergent density estimators required in the process of estimating the desired divergences. In fact, only the empirical CDF is used and the density model being used in the estimator is completely based on the slopes of the piecewise linear interpolation of the empirical CDF. This empirical density model is far from being convergent as we can see from the following figures in Figure 2, which shows the calculated slopes (in blue) for $N=10,100,1000,\mathrm{10,000}$ data samples from a normal distribution against the ground-truth normal densities (in red), plotted in log scale. Clearly, the empirical density model does not converge to the true densities.

Perez-Cruz [1] also provided an empirical KL estimator for multivariate distribution samples. The estimator is based on a nearest-neighbor approach. For each sample $x_i$ in $\mathcal{X}$, where the dimension of the sample is d, let:

${\widehat{p}}_k\left(x_i\right)=\frac{k}{n-1}\frac{\Gamma (d/2+1)}{{\pi }^{d/2}{r}_k{\left(x_i\right)}^d},{\widehat{q}}_k\left(x_i\right)=\frac{k}{m}\frac{\Gamma (d/2+1)}{{\pi }^{d/2}{s}_k{\left(x_i\right)}^d},$

where $r_k\left(x_i\right),s_k\left(x_i\right)$ are, respectively, the Euclidean distance to the k-th nearest neighbor of $x_i$ in $\mathcal{X}\backslash x_i$ and ${\mathcal{X}}^{\prime }$, and $\frac{{\pi }^{d/2}}{\Gamma (d/2+1)}$ is the volume of the unit ball in $R^d$. Ref. [1] continued to show that the random variable $\frac{p\left(x\right)}{{\widehat{p}}_k\left(x\right)}$ converges to an independent Gamma(k,k) random variable which has mean 1 and variance $\frac{1}{k}$ for each selected $k=1,2,3,\dots$, where x is sampled from P. Therefore, they proposed the following estimator:

${\widehat{D}}_k\left(P\right|\left|Q\right)=\frac{1}{n}\sum _{i=1}^{n}log\frac{{\widehat{p}}_k\left(x_i\right)}{{\widehat{q}}_k\left(x_i\right)}=\frac{d}{n}\sum _{i=1}^{n}log\frac{r_k\left(x_i\right)}{s_k\left(x_i\right)}+log\frac{m}{n-1}.$

It was shown that, since $\frac{1}{n}{\sum }_{i=1}^{n}log\frac{p\left(x_i\right)}{{\widehat{p}}_k\left(x_i\right)}$ and, consequently, $\frac{1}{n}{\sum }_{i=1}^{n}log\frac{q\left(x_i\right)}{{\widehat{q}}_k\left(x_i\right)}$ converges to:

$\frac{1}{(k-1)!}{\int }_0^{\infty }{\left(kx\right)}^{k-1}logx{e}^{-kx}kdx=\frac{1}{(k-1)!}{\int }_0^{\infty }{z}^{k-1}logz{e}^{-z}dz-logk,$

then ${\widehat{D}}_k\left(P\right|\left|Q\right)\to D\left(P\right|\left|Q\right)$ almost surely.

4. Empirical Squared Hellinger Distance Estimator of Continuous Distributions

4.1. Estimator for 1D Data

Following Perez-Cruz [1], we have defined a similar estimator for Hellinger affinity using empirical CDFs. Let $\mathcal{X}={\left\{x_i\right\}}_{i=1}^n,{\mathcal{X}}^{\prime }={\left\{x_j^{\prime }\right\}}_{j=1}^m$ be iid samples from P and Q in increasing order. Recall that the definition of the empirical CDFs of P and Q are, respectively,

$P_e\left(x\right)=\frac{1}{n}\sum _{i=1}^{n}U(x-x_i),$

where $U\left(x\right)$ is a unit-step function with $U\left(0\right)=0.5$. The continuous piece-wise linear interpolation of the empirical CDF of P (and Q) is denoted as $P_c\left(x\right)$ (and $Q_c\left(x\right)$, respectively). It is zero for anything smaller than a joint lower bound $x_0<inf\{\mathcal{X},{\mathcal{X}}^{\prime }\}$ of the data samples from $P,Q$, and is one for anything greater than or equal to a joint upper bound $x_{n+1}>sup\{\mathcal{X},{\mathcal{X}}^{\prime }\}$ of the data samples from $P,Q$; everywhere in the middle, it is defined as:

$P_c\left(x\right)=a_{i}x+b_i,x_{i-1}<x<x_i,$

where coefficients $a_i,b_i$ are set so that $P_c\left(x\right)$ matches the values of $P_e\left(x\right)$ at the sampled values $x_i,i=1,\dots ,n$. $Q_c\left(x\right)$ is defined similarly. These empirical CDFs converge uniformly and are independent of the distribution of their CDFs.

Our estimator for the squared Hellinger distance is based on estimating the Hellinger affinity, which is directly related to the quantity of interest by:

$A(P,Q)=1-H^2(P,Q)={\int }_x\sqrt{p\left(x\right)q\left(x\right)}dx.$

The new estimator for Hellinger affinity is

$\widehat{A}(P,Q)=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\delta Q_c\left(x_i\right)}{\delta P_c\left(x_i\right)}},$

where $\delta P_c\left(x_i\right)=(P_c\left(x_i\right)-P_c(x_i-ϵ))/ϵ$ for any $ϵ<{min}_i\{x_i-x_{i-1}\}$ denotes the left slope of $P_c$ at $x_i$ and, similarly, $\delta Q_c\left(x_i\right)$ denotes the left slope of $Q_c$ at $x_i$.

We next claim and prove that $\widehat{A}$ converges to a scalar multiple of the true Hellinger affinity $\widehat{A}\to \frac{\pi }{4}A$. To justify the use of this bias correction constant we need to prove that it results from terms we get from rewriting the estimator:

$\widehat{A}(P,Q)=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q_c\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P_c\left(x_i\right)/\Delta x_i}}$

$=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P\left(x_i\right)/\Delta x_i}}\sqrt{\frac{\Delta Q_c\left(x_{mi}^{\prime }\right)}{\Delta Q\left(x_{mi}^{\prime }\right)}}\sqrt{\frac{\Delta P\left(x_i\right)}{\Delta P_c\left(x_i\right)}}=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P\left(x_i\right)/\Delta x_i}}\frac{\sqrt{n\Delta P\left(x_i\right)}}{\sqrt{m\Delta Q\left(x_{mi}^{\prime }\right)}},$

where $\Delta x_i=x_i-x_{i-1},$ $\Delta P_c\left(x_i\right)=P_c\left(x_i\right)-P_c\left(x_{i-1}\right),$ $\Delta P\left(x_i\right)=P\left(x_i\right)-P\left(x_{i-1}\right),$ $\Delta x_{mi}^{\prime }=min\left\{x_j^{\prime }\right|x_j^{\prime }\ge x_i\}-max\left\{x_j^{\prime }\right|x_j^{\prime }<x_i\},$ $\Delta Q_c\left(x_{mi}^{\prime }\right)=Q_c(min\left\{x_j^{\prime }\right|x_j^{\prime }\ge x_i\})-Q_c(max\left\{x_j^{\prime }\right|x_j^{\prime }<x_i\left\}\right)$ and $\Delta Q\left(x_{mi}^{\prime }\right)=Q(min\left\{x_j^{\prime }\right|x_j^{\prime }\ge x_i\})-Q(max\left\{x_j^{\prime }\right|x_j^{\prime }<x_i\left\}\right)$.

Notice that the first (square root) term in the sum converges almost surely to $\sqrt{\frac{q\left(x_i\right)}{p\left(x_i\right)}}$. We need to show that the above empirical sum converges almost surely to $C{\int }_x\sqrt{p\left(x\right)q\left(x\right)}dx$, where the constant $C=\frac{\pi }{4}$ is derived from the second term, using similar arguments as Perez-Cruz [1] through waiting time distributions between two consecutive samples from a uniform distribution between 0 and 1.

We outline the proof for the constant term below. Similar to Perez-Cruz [1], we know that, given ${\left\{x_i\right\}}_{i=1}^n\sim P$, $n\Delta P\left(x_i\right)\sim Exp\left(1\right)$ and is independent of P (similarly for Q). With this argument, the last expression for $\widehat{A}(P,Q)$ can be rewritten as (where $z_i=n\Delta P\left(x_i\right)$)

$\widehat{A}(P,Q)=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P\left(x_i\right)/\Delta x_i}}\frac{\sqrt{z_i}}{\sqrt{m\Delta Q\left(x_{mi}^{\prime }\right)}}$

$\stackrel{a.s.}{\to }\left(\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P\left(x_i\right)/\Delta x_i}}\sqrt{z_i}\right)\ast \left(\frac{1}{n}\sum _{i=1}^n\frac{1}{\sqrt{m\Delta Q\left(x_{mi}^{\prime }\right)}}\right).$

The first sum converges almost surely to:

$\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{\Delta Q\left(x_{mi}^{\prime }\right)/\Delta x_{mi}^{\prime }}{\Delta P\left(x_i\right)/\Delta x_i}}\sqrt{z_i}\stackrel{a.s.}{\to }{\int }_x{\int }_{z=0}^{\infty }\sqrt{\frac{q\left(x\right)}{p\left(x\right)}}\sqrt{z}{e}^{-z}p\left(x\right)dzdx=\frac{\sqrt{\pi }}{2}A(P,Q).$

The second sum can be rewritten as:

$\frac{1}{n}\sum _{i=1}^n\frac{1}{\sqrt{m\Delta Q\left(x_{mi}^{\prime }\right)}}=\frac{1}{n}\sum _{j=1}^m\frac{n\Delta P_e\left(x_j^{\prime }\right)}{\sqrt{m\Delta Q\left(x_j^{\prime }\right)}}=\frac{1}{m}\sum _{j=1}^m\frac{\Delta P_e\left(x_j^{\prime }\right)/\Delta x_j^{\prime }}{\Delta Q\left(x_j^{\prime }\right)/\Delta x_j^{\prime }}\frac{m\Delta Q\left(x_j^{\prime }\right)}{\sqrt{m\Delta Q\left(x_j^{\prime }\right)}}.$

The last expression converges almost surely to:

$\frac{1}{m}\sum _{j=1}^m\frac{\Delta P_e\left(x_j^{\prime }\right)/\Delta x_j^{\prime }}{\Delta Q\left(x_j^{\prime }\right)/\Delta x_j^{\prime }}\frac{m\Delta Q\left(x_j^{\prime }\right)}{\sqrt{m\Delta Q\left(x_j^{\prime }\right)}}\stackrel{a.s.}{\to }{\int }_x{\int }_{z=0}^{\infty }\frac{p_e\left(x\right)}{q\left(x\right)}\sqrt{z}{e}^{-z}q\left(x\right)dzdx=\frac{\sqrt{\pi }}{2}{\int }_x{p}_e\left(x\right)dx=\frac{\sqrt{\pi }}{2}.$

Notice here that $p_e\left(x\right)$ is a density model but does not need to converge to $p\left(x\right)$ for the above expression to converge to the desired constant.

Combining all previous results, we have shown that $\widehat{A}(P,Q)$ converges almost surely to:

$\widehat{A}(P,Q)\stackrel{a.s.}{\to }\frac{\sqrt{\pi }}{2}\frac{\sqrt{\pi }}{2}A(P,Q)=\frac{\pi }{4}A(P,Q).$

Hence, we obtained the desired constant $C=\frac{\pi }{4}\approx 0.785$. The final estimator for squared Hellinger distance is ${\widehat{H}}^2(P,Q)=1-\frac{4\widehat{A}(P,Q)}{\pi }$.

Notice that Hellinger distance is a symmetric distance metric for any distributions P and Q, hence the estimator above is only one side of the story. Following exactly the same arguments, we can show that the opposite direction estimator,

$\widehat{A}(Q,P)=\frac{1}{m}\sum _{j=1}^m\sqrt{\frac{\delta P_c\left(x_j^{\prime }\right)}{\delta Q_c\left(x_j^{\prime }\right)}},$

also converges almost surely to $\frac{\pi }{4}A(Q,P)$, and since $A(P,Q)=A(Q,P)$ we can obtain a symmetric estimator of Hellinger affinity that converges almost surely to $\frac{\pi }{4}A(Q,P)$:

${\widehat{A}}_S(P,Q)=\frac{\widehat{A}(P,Q)+\widehat{A}(Q,P)}{2}.$

Therefore, we can construct a corresponding estimator for the squared Hellinger distance as:

${\widehat{H}}_S^2(P,Q)=1-\frac{4{\widehat{A}}_S(P,Q)}{\pi },$

which enjoys all of the properties shown above for the two estimators separately. Since the symmetric version uses more information from the two samples, it is supposed to be able to provide better estimates than the two single-sided estimators in terms of the rate of convergence.

4.2. Numerical Experiments

We show asymptotic convergence of the new estimator ${\widehat{H}}^2(P,Q)=1-\frac{4\widehat{A}(P,Q)}{\pi }$, and its symmetric version, to the true $H^2$ value as the data sample size grows in the below experiments. In each of the experiments, we took two distributions of the same family and compared the estimated squared Hellinger distance value against the ground truth value. We plotted mean estimated values for sample size $N=M=10,32,100,316,1000,3162,\mathrm{10,000}$ (x-axis) used for each pair of distributions over 100 instances, and we also plotted the $95\%$ confidence interval of the estimates. For each experiment, the squared Hellinger distance estimators ${\widehat{H}}^2(P,Q),{\widehat{H}}^2(Q,P)$ are plotted in red and blue, and the symmetric squared Hellinger distance estimator ${\widehat{H}}_S^2(P,Q)$ is plotted in purple. We also recall the fact that when $P,Q$ are taken to be normal distributions $N({\mu }_1,{\sigma }_1^2),N({\mu }_2,{\sigma }_2^2)$, the squared Hellinger distance has an analytic form:

$H^2(P,Q)=1-\sqrt{\frac{2{\sigma }_1{\sigma }_2}{{\sigma }_1^2+{\sigma }_2^2}}{e}^{-\frac{1}{4}\frac{{({\mu }_1-{\mu }_2)}^2}{{\sigma }_1^2+{\sigma }_2^2}}.$

In the first experiment (Figure 3), $P,Q$ are taken to be normal distributions $N(0,4)$ and $N(1,1)$, respectively. In the second experiment $P,Q$ are taken to be normal distributions $N(0,1)$ and $N(2,1)$, respectively.

In the third experiment (Figure 4), $P,Q$ are taken to be normal distributions $N(0,1)$ and $N(0.01,1)$, respectively. In the fourth experiment, $P,Q$ are taken to be exponential distributions $Exp\left(1\right)$ and $Exp\left(2\right)$, respectively.

In the fifth experiment (Figure 5), $P,Q$ are taken to be uniform distributions $U(0,1)$ and $U(0,2)$, respectively. In the sixth experiment, $P,Q$ are taken to be uniform distributions $U(0,1)$ and $U(0.5,1.5)$, respectively. Notice that the squared Hellinger distance is well-defined for distributions of different support.

In the last two experiments (Figure 6), we considered two distributions from different distribution families. Here, $P=Cauchy(0,1)$ is the standard Cauchy distribution. In the seventh experiment, $Q=N(1,1)$ and in the last experiment $Q=N(0,1)$. The true squared Hellinger distances are computed using numerical integration.

We can observe from the previous experiments that, depending on the distributions, either the estimator ${\widehat{H}}^2(P,Q)$ or the reverse direction estimator ${\widehat{H}}^2(Q,P)$ can turn out to be better, which is a consequence of our choice to take the left slope of the empirical CDF so the relative location of the two distributions will determine which estimator is more accurate. The symmetric squared Hellinger estimator provides a middle ground between the two one-sided estimators and it also exhibits smaller variances.

As mentioned before, the proposed estimator does not use the information of the underlying distribution and does not need to estimate the density first before estimating the squared Hellinger distance. As a comparison with an estimator that knows the distribution, we performed experiments with Gaussian distributions where we could use the sample mean and sample variance to estimate the distributions and then compute the squared Hellinger distance analytically using the estimated parameters. The estimator is constructed as follows,

${\widehat{H}}_{naive}^2(P,Q)=1-\sqrt{\frac{2{\widehat{\sigma }}_1\widehat{{\sigma }_2}}{{\widehat{\sigma }}_1^2+{\widehat{\sigma }}_2^2}}{e}^{-\frac{1}{4}\frac{{({\widehat{\mu }}_1-{\widehat{\mu }}_2)}^2}{{\widehat{\sigma }}_1^2+{\widehat{\sigma }}_2^2}},$

where ${\widehat{\mu }}_1,{\widehat{\sigma }}_1,{\widehat{\mu }}_2,{\widehat{\sigma }}_2$ are sample estimates of mean and standard deviation from the two datasets. This estimator knows extra information about the data coming from Gaussian distributions.

However, as we can see from the plots in Figure 7, the proposed squared Hellinger distance estimator performs similarly to the estimator that knows the distribution family. In the first experiment, the two distributions are $N(0,4),N(2,4)$. In the second experiment, the two distributions are $N(0,1),N(2,1)$. For both plots, we plotted the proposed symmetric squared Hellinger distance estimator in red and the naive estimator using sampled parameters in blue. The upper bound and lower bound of each estimator, performed over 100 iterations, are plotted in dashed lines.

Finally, we consider the setting in Test 8, where P is a standard Cauchy distribution and Q is a standard normal distribution, and we compare the behavior of the empirical squared Hellinger estimator ${\widehat{H}}^2(P,Q)$ with the empirical KL divergence estimator $\widehat{D}\left(P\right|\left|Q\right)$ as in [1]. As mentioned in the discussion in Section 2, for this case, the KL divergence diverges to infinity while the squared Hellinger distance is bounded. With the same experiment setup, we plotted the resulting divergence estimates and confidence intervals for both KL divergence and squared Hellinger distance in Figure 8, where the ground truth $H^2(P,Q)$ value (approximated by numerical integration) is plotted in black and the ground truth $D_{KL}\left(P\right|\left|Q\right)$ value is infinity.

As we can observe from Figure 8, while the empirical squared Hellinger estimator converges to the ground truth value quickly, the empirical KL divergence estimator cannot converge to some value due to the fact that the ground truth value is infinity. This justifies the desirability of considering the squared Hellinger distance, which is always bounded.

4.3. Estimator for Vectorial Data

Utilizing the results proved for the vectorial data case in [1], we propose the following estimator for squared Hellinger distance in multivariate cases (for a chosen k). Similar to the definitions in [1], let the kNN density estimator be defined as:

${\widehat{p}}_k\left(x_i\right)=\frac{k}{n-1}\frac{\Gamma (d/2+1)}{{\pi }^{d/2}{r}_k{\left(x_i\right)}^d},{\widehat{q}}_k\left(x_i\right)=\frac{k}{m}\frac{\Gamma (d/2+1)}{{\pi }^{d/2}{s}_k{\left(x_i\right)}^d},$

where $r_k\left(x_i\right),s_k\left(x_i\right)$ are, respectively, the Euclidean distance to the k-th nearest neighbor of $x_i$ in $\mathcal{X}\backslash x_i$ and ${\mathcal{X}}^{\prime }$. Let:

${\widehat{A}}_k(P,Q)=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{{\widehat{q}}_k\left(x_i\right)}{{\widehat{p}}_k\left(x_i\right)}}=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{(n-1)r_k{\left(x_i\right)}^d}{m{s}_k{\left(x_i\right)}^d}}$

${\widehat{A}}_k(P,Q)=\frac{1}{n}\sum _{i=1}^n\sqrt{\frac{q\left(x_i\right)}{p\left(x_i\right)}}\sqrt{\frac{{\widehat{q}}_k\left(x_i\right)}{q\left(x_i\right)}}\sqrt{\frac{p\left(x_i\right)}{{\widehat{p}}_k\left(x_i\right)}};$

since $\frac{p\left(x\right)}{{\widehat{p}}_k\left(x\right)},\frac{q\left(x\right)}{{\widehat{q}}_k\left(x\right)}$ are independent Gamma(k,k) random variables that are also independent from $P,Q$, we conclude that ${\widehat{A}}_k(P,Q)$ converges almost surely to:

${\widehat{A}}_k(P,Q)\to A(P,Q)\sqrt{k}{\int }_0^{\infty }{z}^{-1/2}\frac{z^{k-1}{e}^{-z}}{(k-1)!}dz\frac{1}{\sqrt{k}}{\int }_0^{\infty }{z}^{1/2}\frac{z^{k-1}{e}^{-z}}{(k-1)!}dz$

$=\frac{\Gamma (k-\frac{1}{2})}{(k-1)!}\frac{\Gamma (k+\frac{1}{2})}{(k-1)!}A(P,Q).$

So, ${\widehat{A}}_k(P,Q)$ converges almost surely to the true Hellinger affinity up to a constant multiplier, similar to the 1D case. Therefore, we propose the following estimator for squared Hellinger distance, which converges almost surely to the true squared Hellinger distance:

${\widehat{H}}_k^2(P,Q)=1-\frac{{\widehat{A}}_k(P,Q)(k-1)!(k-1)!}{\Gamma (k-\frac{1}{2})\Gamma (k+\frac{1}{2})}\to H^2(P,Q).$

Similar to the 1D case, we can extend this estimator to a symmetric version that also shares the desired convergence properties:

${\widehat{H}}_{k,S}^2(P,Q)=\frac{{\widehat{H}}_k^2(P,Q)+{\widehat{H}}_k^2(Q,P)}{2}.$

4.4. Numerical Experiments for Vectorial Data

Similar to the experiment setting in Section 4.2, we show the convergence of the proposed estimators in Section 4.3. Sample size $N=M=10,32,100,316,1000,3162$ (plotted on the x-axis) is used for each pair of distributions. The analytical formula for the squared Hellinger distance for two multivariate Gaussians $N({\mu }_1,{\Sigma }_1),N({\mu }_2,{\Sigma }_2)$ is:

$H^2(P,Q)=1-\frac{|{\Sigma }_1{|}^{1/4}{\left|{\Sigma }_2\right|}^{1/4}}{|\frac{1}{2}{\Sigma }_1+\frac{1}{2}{\Sigma }_2{|}^{1/2}}{e}^{-\frac{1}{8}{({\mu }_1-{\mu }_2)}^T{(\frac{1}{2}{\Sigma }_1+\frac{1}{2}{\Sigma }_2)}^{-1}({\mu }_1-{\mu }_2)}.$

In the experiment in Figure 9, we picked 2D normal distributions P and Q with ${\mu }_1={(0,0)}^T,{\mu }_2={(1,1)}^T,{\Sigma }_1={\Sigma }_2=I_2$. For the proposed k-nearest neighbor estimator, we picked $k=5$. The performance of the proposed estimators and a comparison with the naive estimator are plotted below. The naive estimator estimates the mean and covariance based on the data samples and estimates the squared Hellinger distance based on the analytic formula:

${\widehat{H}}_{naive}^2(P,Q)=1-\frac{|{\widehat{\Sigma }}_1{|}^{1/4}{\left|{\widehat{\Sigma }}_2\right|}^{1/4}}{|\frac{1}{2}{\widehat{\Sigma }}_1+\frac{1}{2}{\widehat{\Sigma }}_2{|}^{1/2}}{e}^{-\frac{1}{8}{({\widehat{\mu }}_1-{\widehat{\mu }}_2)}^T{(\frac{1}{2}{\widehat{\Sigma }}_1+\frac{1}{2}{\widehat{\Sigma }}_2)}^{-1}({\widehat{\mu }}_1-{\widehat{\mu }}_2)}.$

From these results we can observe that, similar to the 1D cases, the symmetric estimator seems to perform the best and is comparable to the naive estimator in terms of convergence.

In general, a larger k leads to a smaller variance in the proposed estimator for multivariate data. To balance the convergence rate with computational cost, we can select k to be around 4 to 6 which converges faster than a smaller k and is also easy to compute. This behavior is shown in Figure 10, where we compared the performance of the proposed estimator using $k=2,3,4,5,6$ for the same experiment setting as above.

Another test we conducted was to check if the squared Hellinger distance estimate behavior in a non-asymptotic sense is similar for two pairs of concentric Gaussians that have the same squared Hellinger distance. For this experiment, we picked the first pair of Gaussians to be $N(0,I)$ and $N(0,4I)$, and the second pair of Gaussians to be $N(0,\frac{1}{2}I)$ and $N(0,2I)$. The squared Hellinger distance between each pair of Gaussians is 0.2 and, since these two pairs correspond to a single coordinate transformation on the sample space, we expect similar behavior of the estimator in terms of convergence on both pairs. The result is shown in Figure 11. As expected, the empirical estimator for vectorial data has very similar convergence behavior for each of the two pairs of Gaussians to the same ground-truth value.

5. Empirical $\mathbf{\alpha }$-Divergence Estimator of Continuous Distributions

5.1. Estimator for 1D Data

We generalized the results obtained before to a family of $\alpha$-divergences to which the squared Hellinger distance belongs. Following Cichoki and Amari [9], we define an $\alpha$-divergence between two probability distributions as:

$D_A^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha (\alpha -1)}{\int }_x\left(p^{\alpha }\left(x\right)q^{1-\alpha }\left(x\right)-\alpha p\left(x\right)+(\alpha -1)q\left(x\right)\right)dx.$

We want to obtain an empirical estimator similar to that in Section 3 that uses only the empirical CDFs of P and Q and estimates this quantity directly for any $\alpha \in (0,1)$. Notice that for $\alpha =0.5$, $D_A^{\alpha }\left(P\right|\left|Q\right)=4H^2(P,Q)$, which corresponds to the squared Hellinger distance.

Notice that we can rewrite the $\alpha$-divergence above as:

$D_A^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha }+\frac{1}{1-\alpha }-\frac{1}{\alpha (1-\alpha )}{\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx.$

Clearly, we are interested in the last quantity, so we only need to have an estimator for that term that converges almost surely.

For this purpose, let us define an estimator:

${\widehat{A}}^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{\delta Q_c\left(x_i\right)}{\delta P_c\left(x_i\right)}\right)}^{1-\alpha }.$

Notice that, in the general cases, the $\alpha$-divergence is not symmetric.

Following similar procedures as in Section 3, we can rewrite the estimator as:

${\widehat{A}}^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{\Delta Q(x_{mi}^{\prime }/\Delta x_{mi}^{\prime })}{\Delta P\left(x_i\right)/\Delta x_i}\right)}^{1-\alpha }{\left(\frac{n\Delta P\left(x_i\right)}{m\Delta Q\left(x_{mi}^{\prime }\right)}\right)}^{1-\alpha }$

$=\frac{1}{n}\sum _{i=1}^n{\left(\frac{\Delta Q(x_{mi}^{\prime }/\Delta x_{mi}^{\prime })}{\Delta P\left(x_i\right)/\Delta x_i}\right)}^{1-\alpha }{\left(\frac{z_i}{m\Delta Q\left(x_{mi}^{\prime }\right)}\right)}^{1-\alpha }.$

This sum converges almost surely to (since the exponential waiting distributions are independent of the data distribution):

$\left(\frac{1}{n}\sum _{i=1}^n{\left(\frac{\Delta Q(x_{mi}^{\prime }/\Delta x_{mi}^{\prime })}{\Delta P\left(x_i\right)/\Delta x_i}\right)}^{1-\alpha }{z}_i^{1-\alpha }\right)(\frac{1}{m}\sum _{j=1}^m{\left(\frac{1}{m\Delta Q\left(x_j^{\prime }\right)}\right)}^{1-\alpha }m\Delta Q\left(x_j^{\prime }\right)\frac{\Delta P_e\left(x_j^{\prime }\right)}{\Delta Q\left(x_j^{\prime }\right)}).$

Following the same arguments as in Section 3, we can show that the proposed estimator converges almost surely to:

${\widehat{A}}^{\alpha }\left(P\right|\left|Q\right)\stackrel{a.s.}{\to }{\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx{\int }_{z=0}^{\infty }{z}^{1-\alpha }{e}^{-z}dz{\int }_{z=0}^{\infty }{z}^{\alpha }{e}^{-z}dz{\int }_x{p}_e\left(x\right)dx$

$=C_{1-\alpha }{C}_{\alpha }{\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx,$

where we define the constants $C_{1-\alpha }={\int }_{z=0}^{\infty }{z}^{1-\alpha }{e}^{-z}dz=\Gamma (2-\alpha ),C_{\alpha }={\int }_{z=0}^{\infty }{z}^{\alpha }{e}^{-z}dz=\Gamma (1+\alpha ),\forall \alpha \in (-1,2)$.

Therefore, we know that the estimator

${\widehat{D}}_A^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha }+\frac{1}{1-\alpha }-\frac{1}{\alpha (1-\alpha )}\frac{{\widehat{A}}^{\alpha }\left(P\right|\left|Q\right)}{C_{\alpha }{C}_{1-\alpha }}$

converges almost surely to the true $\alpha$-divergence value, $D_A^{\alpha }\left(P\right|\left|Q\right)$.

Although the $\alpha$-divergence is not symmetric, it has the property that

$D_A^{\alpha }\left(P\right|\left|Q\right)=D_A^{1-\alpha }\left(Q\right|\left|P\right).$

So, given the same two sample data sets, we can get another estimator for the same quantity based on ${\widehat{D}}_A^{1-\alpha }\left(Q\right|\left|P\right)=\frac{1}{\alpha }+\frac{1}{1-\alpha }-\frac{1}{\alpha (1-\alpha )}\frac{{\widehat{A}}^{1-\alpha }\left(Q\right|\left|P\right)}{C_{\alpha }{C}_{1-\alpha }}$, where we are estimating based on the sampling distribution from Q instead of P. Since ${\widehat{D}}_A^{\alpha }\left(P\right|\left|Q\right),{\widehat{D}}_A^{1-\alpha }\left(Q\right|\left|P\right)$ converges to the same divergence value, we can again create a symmetric estimator based on averaging these two estimators ${\widehat{D}}_{A,S}^{\alpha }\left(P\right|\left|Q\right)=\frac{{\widehat{D}}_A^{\alpha }\left(P\right|\left|Q\right)+{\widehat{D}}_A^{1-\alpha }\left(Q\right|\left|P\right)}{2}$ and it is expected to perform similarly if not better. Lastly, notice that, when $\alpha =0.5$, we obtain $C_{\alpha }=C_{1-\alpha }=\frac{\sqrt{\pi }}{2}$ and ${\widehat{D}}_A^{0.5}\left(P\right|\left|Q\right)=4(1-\frac{4}{\pi }{\widehat{A}}^{0.5}\left(P\right|\left|Q\right))$, which corresponds to the squared Hellinger estimator we have seen in Section 3, scaled by 4.

5.2. Numerical Experiments

We show asymptotic convergence of the new estimator ${\widehat{D}}_A^{\alpha }\left(P\right|\left|Q\right)$, and its symmetric version, to the true $\alpha$-divergence value as the data sample size grows in the below experiments. In each of the below experiments, we took two distributions of the same family and compared the estimated $\alpha$-divergence value against the ground truth value. Mean estimated values for sample size $N=M=10,32,100,316,1000,3162,\mathrm{10,000},\mathrm{31,623}$ (plotted on the x-axis) used for each pair of distributions over 100 instances and we also plotted the $95\%$ confidence interval of the estimates. For each experiment, the $\alpha$-divergence estimators ${\widehat{D}}_A^{\alpha }\left(P\right|\left|Q\right),{\widehat{D}}_A^{1-\alpha }\left(Q\right|\left|P\right)$ are plotted in red and blue, and the symmetric $\alpha$-divergence estimator ${\widehat{D}}_{A,S}^{\alpha }\left(P\right|\left|Q\right)$ is plotted in purple.

In the first experiment, $P,Q$ are taken to be normal distributions $N(0,4)$ and $N(1,1)$, respectively, and $\alpha =0.6$. In the second experiment, $P,Q$ are taken to be normal distributions $N(0,1)$ and $N(2,1)$, respectively, and $\alpha =0.4$. The results are plotted in Figure 12. Notice that for two normal distributions $P\sim N({\mu }_1,{\sigma }_1^2),Q\sim N({\mu }_2,{\sigma }_2^2)$, we have an analytical formula for the $\alpha$-divergence:

$D_A^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha (1-\alpha )}\left(1-\frac{{\sigma }_2^{\alpha }{\sigma }_1^{1-\alpha }}{\sqrt{\alpha {\sigma }_2^2+(1-\alpha ){\sigma }_1^2}}{e}^{-\frac{\alpha (1-\alpha )}{\alpha {\sigma }_2^2+(1-\alpha ){\sigma }_1^2}\frac{{({\mu }_1-{\mu }_2)}^2}{2}}\right).$

Again, we provide a comparison with an estimator that knows the distribution family. We performed experiments with Gaussian distributions where we could use the sample mean and sample variance to estimate the distributions and then compute the $\alpha$-divergences analytically using the estimated parameters. The estimator is constructed as follows:

${\widehat{D}}_{A,naive}^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha (1-\alpha )}\left(1-\frac{{\widehat{\sigma }}_2^{\alpha }{\widehat{\sigma }}_1^{1-\alpha }}{\sqrt{\alpha {\widehat{\sigma }}_2^2+(1-\alpha ){\widehat{\sigma }}_1^2}}{e}^{-\frac{\alpha (1-\alpha )}{\alpha {\widehat{\sigma }}_2^2+(1-\alpha ){\widehat{\sigma }}_1^2}\frac{{({\widehat{\mu }}_1-{\widehat{\mu }}_2)}^2}{2}}\right)$

where ${\widehat{\mu }}_1,{\widehat{\sigma }}_1,{\widehat{\mu }}_2,{\widehat{\sigma }}_2$ are sample estimates of mean and standard deviation from the two datasets. This estimator knows extra information about the data coming from Gaussian distributions. However, as we can see from the plots in Figure 13, the proposed $\alpha$-divergence estimator performs similarly to the estimator that knows the distribution family.

In the first experiment, the two distributions are $N(0,4),N(1,1)$ and $\alpha =0.6$. In the second experiment, the two distributions are $N(0,1),N(2,1)$ and $\alpha =0.4$. For both plots, we plotted the proposed symmetric $\alpha$-divergence estimator in red and the naive estimator using sampled parameters in blue. The upper bound and lower bound of each estimator, performed over 100 iterations, are plotted in dashed lines.

5.3. Estimator for Vectorial Data

Similarly to Section 4.3, we propose $\alpha$-divergence estimators for samples from multivariate distributions. For this purpose, let us define:

${\widehat{A}}_k^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{n}\sum _{i=1}^n{\left(\frac{{\widehat{q}}_k\left(x_i\right)}{{\widehat{p}}_k\left(x_i\right)}\right)}^{1-\alpha }.$

Using similar arguments, we can show that this estimator converges almost surely to:

${\widehat{A}}_k^{\alpha }\left(P\right|\left|Q\right)\to \left(k^{1-\alpha }{\int }_0^{\infty }{z}^{\alpha -1}\frac{z^{k-1}{e}^{-z}}{(k-1)!}dz\right)\left(k^{\alpha -1}{\int }_0^{\infty }{z}^{1-\alpha }\frac{z^{k-1}{e}^{-z}}{(k-1)!}dz\right){\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx$

$=\frac{\Gamma (k+\alpha -1)}{(k-1)!}\frac{\Gamma (k-\alpha +1)}{(k-1)!}{\int }_x{\left(\frac{q\left(x\right)}{p\left(x\right)}\right)}^{1-\alpha }p\left(x\right)dx.$

Therefore, we propose the following estimator for $\alpha$-divergences, which converges almost surely:

${\widehat{D}}_{A,k}^{\alpha }\left(P\right|\left|Q\right)=\frac{1}{\alpha }+\frac{1}{1-\alpha }-\frac{1}{\alpha (1-\alpha )}\frac{{\widehat{A}}_k^{\alpha }\left(P\right|\left|Q\right)(k-1)!(k-1)!}{\Gamma (k+\alpha -1)\Gamma (k-\alpha +1)}\to D_A^{\alpha }\left(P\right|\left|Q\right).$

Similarly, we can extend this estimator to a symmetric version, for any fixed k:

${\widehat{D}}_{A,k,S}^{\alpha }\left(P\right|\left|Q\right)=\frac{{\widehat{D}}_{A,k}^{\alpha }\left(P\right|\left|Q\right)+{\widehat{D}}_{A,k}^{1-\alpha }\left(Q\right|\left|P\right)}{2}.$

As a remark, for the vectorial case, the above kNN density-based empirical estimator for $\alpha$-divergences (and the squared Hellinger distance in Section 4.3 as a special case) agree with the estimators proposed in [13], although the proof of convergence differs. Nonetheless, the univariate estimators we proposed in Section 4.1 and Section 5.1 are different from trivial reductions of the kNN-based estimators in Section 4.3 and Section 5.3 when taking $d=1$ and $k=1$.

6. Limitation of the Proposed Methodologies and Uniqueness of the $\mathbf{\alpha }$-Divergences

6.1. Failure of a Similar Estimator for Total Variation Distance

As we have shown so far, by using the trick of waiting time distributions, we can bias-correct an empirical mean type estimator to produce an almost-sure convergence estimator for KL divergence, squared Hellinger distance, and in general the $\alpha$-divergences. However, the same kind of trick does not work for other f-divergences that have an f-function without certain desired properties such as $f\left(ab\right)=f\left(a\right)+f\left(b\right)$ for KL divergence or $f\left(ab\right)=f\left(a\right)f\left(b\right)$ for Hellinger affinity, which we shall discuss in more detail later. As a simple demonstration, consider the Total Variation Distance (TVD), which, for two continuous distributions P and Q, is defined as:

$TVD(P,Q)=\frac{1}{2}{\int }_x|p\left(x\right)-q\left(x\right)|dx.$

Notice that the TVD is always bounded between 0 and 1.

We considered paired distributions in two different families in 1D, namely normal distributions and exponential distributions. For different choices of parameters, we plotted the performance of a biased estimator using the empirical CDFs against the true TVD value. For every parameter setting, we looked at the case where $N=M=\mathrm{10,000}$ and averaged over 100 instances. The estimator is defined as:

$\widehat{TVD}(P,Q)=\frac{1}{n}\sum _{i=1}^n\frac{1}{2}|\frac{\delta Q_c\left(x_i\right)}{\delta P_c\left(x_i\right)}-1|$

Specifically, for the normal distributions, we fixed ${\mu }_1=0,{\sigma }_1=1,{\sigma }_2=1$ and varied ${\mu }_2$ from 0 to 5. For the exponential distributions, we fixed ${\lambda }_1=0.1$ and varied ${\lambda }_2$ from 0.1 to 7. This generated a range of true TVD values that are spaced between 0 and 1 for each distribution family. Figure 14 plots the biased estimator values (on the y-axis) against true TVD values (on the x-axis) for pairs of normal distributions $P,Q$ in blue and pairs of exponential distributions $P,Q$ in red. The confidence intervals are also plotted. We observe that, for the same true TVD values, the biased estimator produced different values for different distribution families, where the relationship looks nonlinear and depends on the distribution family itself. This is an indication that the proposed estimator cannot be uniformly corrected with a simple additive and/or multiplicative constant as we performed for squared Hellinger distance (and in general $\alpha$-divergences) and [1] for KL divergences. Therefore, we conclude that, so far, the proposed methodologies work for KL divergence, squared Hellinger distance, and in general $\alpha$-divergences only, but cannot be extended to the general f-divergences in a straightforward way.

6.2. Uniqueness of $\alpha$-Divergences

We provide a more detailed explanation as to why the $\alpha$-divergences are the unique family of f-divergences that can be estimated using our type of estimator based on waiting time random variable transformations. Take the vectorial case for example, where we construct kNN empirical density estimates for the probability densities ${\widehat{p}}_k,{\widehat{q}}_k$; for an estimator that is based on these estimates to work for an f-divergence, we would require the f-divergence to be computable through an affinity term as an integration of the form ${\int }_{x}f\left(\frac{p\left(x\right)}{q\left(x\right)}\right)p\left(x\right)dx$ or ${\int }_{x}f\left(\frac{p\left(x\right)}{q\left(x\right)}\right)q\left(x\right)dx$ up to some constant terms, and we require that the affinity generating functions f satisfy a functional form that can be separated as either $f\left(ab\right)=g\left(a\right)+h\left(b\right)$ or $f\left(ab\right)=g\left(a\right)h\left(b\right)$ for some functions g and h. This restriction is made because, as we have seen for KL or $\alpha$-divergence estimators, we rely on the independence property of the waiting time random variables, hence we can separate the empirical sums into three terms which converge separately and show the estimator to converge asymptotically up to additive or multiplicative bias constants. Let us examine these two types of restrictions on f.

For f to satisfy $f\left(xy\right)=g\left(x\right)+h\left(y\right),\forall x,y>0$, we can see that f is equivalent to g and h in the sense that they differ by a constant. Differentiating the previous equation with respect to x and setting $x=1$ we would get:

$f^{\prime }\left(y\right)=\frac{c}{y},$

where $c=g^{\prime }\left(1\right)$ is a constant. The unique family of solutions to this condition is $f\left(x\right)=clogx$ up to some additive constants. This obviously corresponds to KL divergence and reverse KL divergence when integrated against P and Q, respectively.

For the other case where $f\left(xy\right)=g\left(x\right)h\left(y\right),\forall x,y>0$, let us consider differentiating both sides with respect to x; this gives:

$f^{\prime }\left(xy\right)=g^{\prime }\left(x\right)\frac{h\left(y\right)}{y}.$

Taking log on both sides and let $l=log{f}^{\prime },m=log{g}^{\prime }$:

$l\left(xy\right)=m\left(x\right)+log\frac{h\left(y\right)}{y}.$

Now take the derivative with respect to x again and set $x=1$, we get:

$l^{\prime }\left(y\right)=\frac{c}{y},$

where $c=m^{\prime }\left(1\right)$ is a constant. The unique family of solutions satisfying the last condition is $l\left(y\right)=clogy+C$ and hence $f\left(y\right)=a{y}^b$ is a general solution up to some additive constant. Without loss of generality, we can see that this corresponds uniquely to the affinity term of interest of the family of alpha divergences where $f\left(y\right)=y^{\alpha },\forall \alpha \in (0,1)$, and up to some constant terms.

Since KL and reverse KL divergence are limits of the $\alpha$-divergences at two endpoints, we can conclude that the unique family of f-divergences that can be estimated based on the proposed estimators using waiting time random variables are the $\alpha$-divergences. There is an interesting connection, pointed out by Amari [14], that states that $\alpha$-divergence is the unique intersection between f-divergences and decomposable Bregman divergences on the space on positive measures. Notice that if restricted to the space of probability measures then the intersection reduces to only the KL divergences. Although the result does not directly connect to the uniqueness of $\alpha$-divergences being estimable through our proposed methodologies, the proof technique that justifies the functional forms of the $\alpha$-divergence being the unique f-divergence that allows a decomposition into the Bregman divergence dual functions up to some nonlinear coordinate transformation is very similar to what we carried out above and reaches the same conclusion—that the function f must take on a power function form that corresponds to an $\alpha$-divergence and at the limit of $\alpha$ becomes logarithm functions that correspond to KL and reverse KL divergences.

7. Applications

The proposed estimator finds interesting applications in statistical estimation theory, clustering algorithms, visualization/embedding algorithms, and possibly online learning algorithms. We next describe a few such examples.

7.1. Bounding the Neyman–Pearson Region by Hellinger Affinity

We show that the Neyman–Pearson region contains one convex region determined by the Hellinger affinity, which is contained in another. These inclusion relations generalize the classical inequalities between total variation and Hellinger distance. Deploying our estimator for Hellinger affinity ${\widehat{A}}_S(P,Q)$, we can approximately bound the Neyman–Pearson region.

Our results (see Appendix A for more details) show that, with two distributions $p,q$, and with $s,t>0$ (which can be chosen so that $s+t=2$ in standard case), the Neyman–Pearson region for type I ($\alpha \left(E\right)$) and type II ($\beta \left(E\right)$) errors satisfies the following relation with the total variation distance for optimal choice of event $E^{★}$:

$s\alpha \left(E^{★}\right)+t\beta \left(E^{★}\right)=\frac{t+s}{2}-\frac{1}{2}\int |sp-tq|d\mu ,$

and can hence be bounded by the following inequalities where $\rho (p,q)$ is the Hellinger affinity:

$\frac{s+t}{2}-\sqrt{{\left(\frac{s+t}{2}\right)}^2-st\rho {(p,q)}^2}\le s\alpha \left(E^{★}\right)+t\beta \left(E^{★}\right)\le \sqrt{st}\rho (p,q).$

Hence, by substituting our symmetric estimator for the Hellinger affinity term ${\widehat{A}}_S(p,q)\approx \rho (p,q)$, we can approximately bound the Neyman–Pearson region given two samples from distributions p and q,

$\frac{s+t}{2}-\sqrt{{\left(\frac{s+t}{2}\right)}^2-st{\widehat{A}}_S{(p,q)}^2}⪅s\alpha \left(E^{★}\right)+t\beta \left(E^{★}\right)⪅\sqrt{st}{\widehat{A}}_S(p,q).$

If we are dealing with multivariate distributions, then the appropriate multivariate Hellinger affinity estimator from Section 4.3 can be used to approximately bound the Neyman–Pearson region. As a remark, we observe that there is no provable general relationship between Kullback–Leibler divergence or the rest of the $\alpha$-divergences (besides Hellinger distance) with the Neyman–Pearson regions.

7.2. Estimating Eigenvalues of the Matrix Pencil for Inference in the Family of Concentric Gaussians

Consider two multivariate distributions from the concentric Gaussian family $P=N(0,C_1^2),Q=N(0,C_2^2)$, where $C_1^2,C_2^2\in R^{d\times d}$. It can be shown that any meaningful statistical inference function on the two covariance matrices should satisfy $\varphi (C_1^2,C_2^2)=\varphi (I,\Lambda )$, where $\Lambda$ is the diagonal matrix with diagonal entries ${\lambda }_1,\dots ,{\lambda }_d$ being the eigenvalues of the matrix $C_1^{-1}{C}_2^2{\left(C_1^{-1}\right)}^{*}$; see Appendix C for more details.

Since $\Lambda$ is diagonal and I is simply the identity matrix, we can write $\varphi (C_1^2,C_2^2)=h({\lambda }_1,\dots ,{\lambda }_d)$. Hence, any inference we can make on the two concentric Gaussians will depend only on sufficient statistics, which are the eigenvalues ${\lambda }_1,\dots ,{\lambda }_d$. In the case of Hellinger affinity (and in general affinities for $\alpha$-divergences), we can write it as $\varphi (C_1^2,C_2^2)=h(1,{\lambda }_1)\times \dots \times h(1,{\lambda }_d)$, where $h(1,{\lambda }_i),\forall i=1,\dots ,d$ is the affinity calculated based on two univariate Gaussian distributions $N(0,1)$ and $N(0,{\lambda }_i)$. For example, we have analytic formulas for the affinity term of the $\alpha$-divergence family between such univariate Gaussian distributions:

$h_{\alpha }(1,\lambda )=\sqrt{\frac{{\lambda }^{\alpha }}{\alpha \lambda +(1-\alpha )}}.$

Then, we have, for the d-dimensional multivariate concentric Gaussians,

$A^{\alpha }\left(P\right|\left|Q\right)={\varphi }_{\alpha }(C_1^2,C_2^2)=\sqrt{\frac{{\left({\prod }_{i=1}^d{\lambda }_i\right)}^{\alpha }}{{\prod }_{i=1}^d(\alpha {\lambda }_i+1-\alpha )}}.$

Now, given d distinct values of ${\alpha }_1,\dots ,{\alpha }_d$, the affinity values $A^{{\alpha }_1}\left(P\right|\left|Q\right),\dots ,A^{{\alpha }_d}\left(P\right|\left|Q\right)$ can be used to determine the eigenvalues ${\lambda }_1,\dots ,{\lambda }_d$. Since our proposed estimator for vectorial $\alpha$-affinities ${\widehat{A}}_k^{\alpha }\left(P\right|\left|Q\right)$ converges up to a multiplicative constant, we can use the estimated values for $\alpha$-affinities corresponding to d different values of $\alpha ={\alpha }_1,\dots ,{\alpha }_d$ to estimate the eigenvalues ${\lambda }_1,\dots ,{\lambda }_d$ by solving a system of d equations. The estimated values ${\widehat{\lambda }}_1,\dots ,{\widehat{\lambda }}_d$ can be then used for any inference problems on these two probability distributions and they are sufficient for inference. This significantly reduces the noise in estimating the entire covariance matrices $C_1^2,C_2^2$ when the data come from high dimensions where we could have an over-parametrization problem.

7.3. Stock Clustering and Visualization

We next describe a simple application to stock segmentation in a portfolio allocation setting. Consider N stocks with T historic dates. Let ${\left\{r_{i,t}\right\}}_{i\in \left[N\right],t\in \left[T\right]}$ denote the returns of each stock on each date. Let ${\left\{R_i\right\}}_{i\in \left[N\right]}$ denote the random variable standing for the returns of each stock, which is composed of data ${\left\{r_{i,t}\right\}}_{t\in [T}$. To cluster this universe of stocks into K distinct groups, we can first use the Hellinger distance estimator ${\widehat{H}}_S(R_i,R_j)$ for a pair of stocks $\forall i\ne j\in \left[N\right]$. Since the estimator is symmetric, we would arrive at a symmetric distance matrix denoted by $D_H$. It is also possible to combine the Hellinger distance with a correlation distance metric through some transformations. After obtaining the distance matrix (or an affinity matrix by subtracting it from 1), we can deploy any desired clustering algorithm on it. The result would be K clusters of stocks that are grouped by similarity in the chosen distance sense. We can also add another step, which is to repair the distance matrix before clustering. There is the possibility that the distance matrix estimated using the proposed estimator does not exactly correspond to a metric, which means some groups of stocks may violate the triangle law in a metric. We can apply a simple sparse metric repair algorithm, see, for example, [15]. The resulting clustering can be helpful for portfolio allocation strategies since we can build sub-strategies inside each cluster and merge them together.

Another example using the same distance matrix constructed from sample data is in visualization algorithms such as FATE [16], which allow for the input of a precomputed distance/affinity matrix specifying the dataset. The visualization algorithm uses the input distance to compute embeddings in lower dimensions that preserve the local/global structures of the dataset and can be useful in many subsequent applications. Here, our estimator can also serve to compute the input distance matrix on sample data from N entities using the Hellinger distance or $\alpha$-divergences as the distance metric. This could also be used in conjunction with a metric repair algorithm to adjust for the biases and errors in empirical estimators.

7.4. Other Applications

Lastly, we suspect that the proposed estimator can find interesting applications in UCB-type algorithms in multi-armed bandit frameworks, where the estimated pairwise Hellinger distances/$\alpha$-divergences for sample distributions from different arms can be used to eliminate arms that fall outside of the confidence region balls around the top arms historically. We leave these open problems as future works.

8. Conclusions

We have proposed an estimator for the Hellinger affinity, and hence the squared Hellinger distance, between samples from two distributions based solely on the empirical CDF without the need to estimate the densities themselves. We have proven its almost-sure convergence to the true squared Hellinger distance and have constructed a symmetric version of this estimator. We showed the convergence behavior using several experiments where we observed that the symmetric estimator constructed from averaging the two one-sided estimators for the squared Hellinger distance turned out to be a favorable choice due to accuracy in general and smaller variances. We then extended the estimator to a family of $\alpha$-divergences, where similar properties hold up to small modifications. For each choice of $\alpha$, we also showed how to construct a symmetric version of the estimator. We also extended respective estimators to work with multivariate data in higher dimensions using k-nearest-neighbor-based estimators. Numerical examples are given to show the convergence of our proposed estimators. We conclude that the $\alpha$-divergence family is the unique f-divergences that can be estimated consistently using the proposed methodologies. Our proposed estimators can be applied to approximately bounding the Neyman–Pearson region of a statistical test, among many other applications.

Author Contributions

Methodology, R.D.; Formal analysis, R.D.; Writing—original draft, R.D.; Writing—review & editing, A.M.; Visualization, R.D.; Supervision, A.M. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Empirical KLD estimator for two normal distributions.

Figure 2. Failure of empirical PDF estimator.

Figure 3. Empirical squared Hellinger estimator tests between 1D normal distributions.

Figure 4. Empirical squared Hellinger estimator tests between 1D distributions.

Figure 5. Empirical squared Hellinger estimator tests between 1D distributions.

Figure 6. Empirical squared Hellinger estimator tests between 1D Cauchy and normal distributions.

Figure 7. Comparisons between empirical and naive estimators for 1D normal distributions.

Figure 8. Comparison of empirical [Forumla omitted. See PDF.] estimator against empirical [Forumla omitted. See PDF.] estimator, [Forumla omitted. See PDF.][Forumla omitted. See PDF.].

Figure 9. Vectorial squared Hellinger estimator ([Forumla omitted. See PDF.]) tests on 2D normal distributions.

Figure 10. Comparison of kNN-based squared Hellinger distance estimators.

Figure 11. Two pairs of concentric Gaussians with invariant squared Hellinger distance.

Figure 12. [Forumla omitted. See PDF.]-divergence estimator tests on 1D normal distributions.

Figure 13. Comparisons between empirical and naive estimators for 1D normal distributions.

Figure 14. Raw empirical TVD estimator.

Appendix A. Shannon Entropy Estimator for 1D and Vectorial Data

Another simple extension of the methodologies in this work provides us with a convergent estimator for the Shannon entropy defined as:$$H\left(P\right)=-{\int }_{x}p\left(x\right)logp\left(x\right)dx$$

Here given 1D data samples ${\left\{x_i\right\}}_{i=1}^n$ from distribution P, we propose the following estimator:$$\widehat{H}\left(P\right)=-\frac{1}{n}\sum _{i=1}^{n}log\delta P_c\left(x_i\right)$$ where $\delta P_c\left(x_i\right)$ are as defined in Section 4.1. It can be shown that:$$\widehat{H}\left(P\right)=-\frac{1}{n}\sum _{i=1}^{n}logp\left(x_i\right)+\frac{1}{n}\sum _{i=1}^{n}log\frac{\Delta P\left(x_i\right)}{\Delta P_c\left(x_i\right)}=\frac{1}{n}\sum _{i=1}^{n}logn\Delta P\left(x_i\right)-\frac{1}{n}\sum _{i=1}^{n}logp\left(x_i\right)$$ This suggests that $\widehat{H}\left(P\right)\to C+H\left(P\right)$ where $C={\int }_0^{\infty }logz{e}^{-z}dz\approx -0.5772$ is the Euler-Mascheroni constant. Hence we conclude that $\widehat{H}\left(P\right)+0.5772$ converges almost surely to $H\left(P\right)$, the true Shannon entropy of P.

Similarly, for vectorial data in d-dimensions, we define the estimator based on the k-nearest neighbor for a fixed k:$${\widehat{H}}_k\left(P\right)=-\frac{1}{n}\sum _{i=1}^{n}log{\widehat{p}}_k\left(x_i\right)$$ where ${\widehat{p}}_k\left(x_i\right)$ is as defined in Section 4.3. By a similar argument, we show that:$${\widehat{H}}_k\left(P\right)=-\frac{1}{n}\sum _{i=1}^{n}logp\left(x_i\right)+\frac{1}{n}\sum _{i=1}^{n}log\frac{p\left(x_i\right)}{{\widehat{p}}_k\left(x_i\right)}\to H\left(P\right)+\frac{1}{(k-1)!}{\int }_0^{\infty }{z}^{k-1}logz{e}^{-z}dz-logk$$ Since the integral ${\int }_0^{\infty }{z}^{k-1}logz{e}^{-z}dz$ evaluates to ${\sum }_{j=1}^{k-1}\frac{1}{j}+C$, we conclude that ${\widehat{H}}_k\left(P\right)+logk-{\sum }_{j=1}^{k-1}\frac{1}{j}+0.5772$ converges almost surely to $H\left(P\right)$.

On a related note, a class of estimators of the Rényi and Tsallis entropies for multidimensional densities has been studied in [17], which is also based on computing the k-nearest neighbor distances from empirical samples.

Appendix B. Hellinger Affinity and Neyman–Pearson Region

Following [18], we define the squared Hellinger distance and Hellinger affinity as:$$H^2(p,q)=\frac{1}{2}\int {|\sqrt{p}-\sqrt{q}|}^{2}d\mu =1-\rho (p,q)$$ Then the Type I and Type II errors for any event E used as a test for distribution by q, are given by $$\alpha \left(E\right)={\int }_{E}pd\mu$$ $$\beta \left(E\right)=1-{\int }_{E}qd\mu$$ and we have the inequality for any non-negative $s,t$:$$s\alpha \left(E\right)+t\beta \left(E\right)\ge \frac{t+s}{2}-\frac{1}{2}\int |sp-tq|d\mu$$ This is because $s\alpha \left(E\right)+t\beta \left(E\right)=t-{\int }_E(tq-sp)d\mu \ge t-{\int }_E|sp-tq|d\mu$, and $s\alpha \left(E\right)+t\beta \left(E\right)=s-{\int }_{E^c}(sp-tq)d\mu \ge s-{\int }_{E^c}|sp-tq|d\mu$, where $E^c$ is the complement of event E. Combining these results gives the aforementioned inequality, which can be seen as a generalization of the classic case when $s=t=1$, see Chapter 13 of [18]. For an optimal $E^{★}(s,t)$, an event for which this holds with equality (for example when $E^{★}(s,t)$ is the support of $sp-tq<0$), we have:$$s\alpha \left(E^{★}(s,t)\right)+t\beta \left(E^{★}(s,t)\right)=\frac{t+s}{2}-\frac{1}{2}\int |sp-tq|d\mu$$ and $(\alpha \left(E^{★}(s,t)\right),\beta \left(E^{★}(s,t)\right))$ is the point on the Neyman–Pearson boundary with supporting line $s\alpha \left(E\right)+t\beta \left(E\right)\ge s\alpha \left(E^{★}(s,t)\right)+t\beta \left(E^{★}(s,t)\right)$. Since the Neyman–Pearson region is convex, the family of f-divergences $\frac{1}{2}\int |sp-tq|d\mu$, where say $s+t=2$, obtains a complete description of the Neyman–Pearson region.

Additionally, we can relate the Neyman–Pearson region to the Hellinger distance by using the Hellinger affinity. It is convenient to compute $$\frac{1}{2}\int {|\sqrt{sp}+\sqrt{tq}|}^{2}d\mu =\frac{s+t}{2}+\sqrt{st}\rho (p,q)$$ $$\frac{1}{2}\int {|\sqrt{sp}-\sqrt{tq}|}^{2}d\mu =\frac{s+t}{2}-\sqrt{st}\rho (p,q)$$ where the first term is from the conservation of probability. Then we use this to simplify the chain of inequalities:$$\frac{1}{2}\int |\sqrt{sp}-\sqrt{tq}{|}^{2}d\mu \le \frac{1}{2}\int |\sqrt{sp}-\sqrt{tq}|(\sqrt{sp}+\sqrt{tq})d\mu =\frac{1}{2}\int |sp-tq|d\mu$$ and $$\frac{1}{2}\int |sp-tq|d\mu =\frac{1}{2}\int |\sqrt{sp}-\sqrt{tq}|(\sqrt{sp}+\sqrt{tq})d\mu$$ $$\le (\frac{1}{2}\int |\sqrt{sp}-\sqrt{tq}{{|}^{2}d\mu )}^{\frac{1}{2}}{(\frac{1}{2}\int {(\sqrt{sp}+\sqrt{tq})}^{2}d\mu )}^{\frac{1}{2}}$$ $$=\sqrt{{\left(\frac{s+t}{2}\right)}^2-st\rho {(p,q)}^2}$$ (where we used the Cauchy-Schwarz inequality) to the chain of inequalities:$$\frac{s+t}{2}-\sqrt{st}\rho (p,q)\le \frac{1}{2}\int |sp-tq|d\mu \le \sqrt{{\left(\frac{s+t}{2}\right)}^2-st\rho {(p,q)}^2}$$ We obtain then for the Neyman–Pearson boundary the upper and lower bounds in terms of the Hellinger affinity:$$\frac{s+t}{2}-\sqrt{{\left(\frac{s+t}{2}\right)}^2-st\rho {(p,q)}^2}\le s\alpha \left(E^{★}\right)+t\beta \left(E^{★}\right)\le \sqrt{st}\rho (p,q)$$ This bound has the gigantic advantage of also bounding the Neyman–Pearson region for joint distributions, such as the result of i.i.d. samples. There is no general relationship between Kullback–Leibler divergence and Neyman–Pearson regions since, say, for the Bernoulli family we can have a sequence of pairs $\{p_k,q_k\}$ such that $H^2(p_k,q_k)\to 0$ but $D_{KL}(p_k\left|\right|q_k)+D_{KL}(q_k\left|\right|p_k)\to \infty$.

Appendix C. Sufficient Information Eigenvalues for Inference between Concentric Gaussians

Let $P=N(0,C_1^2),Q=N(0,C_2^2)$. For any meaningful statistical inference function on $P,Q$ we require $$\varphi (C_1^2,C_2^2)=\varphi (X{C}_1^2{X}^{*},X{C}_2^2{X}^{*})$$ where X is a coordinate transformation. Writing the QR decomposition for $X=QR$ and let R be chosen so that $R{C}_1^2{R}^{*}=I$ (where Q is unitary), which is equivalent to the Cholesky factorization $C_1^{-2}=R^{*}R$. Define $C_1=R^{-1}$ then we have $$\varphi (C_1^2,C_2^2)=\varphi (I,Q\left(C_1^{-1}{C}_2^2{\left(C_1^{-1}\right)}^{*}\right)Q^{*})$$

We choose Q so that it is the eigenvectors of the Hermitian matrix $C_1^{-1}{C}_2^2{\left(C_1^{-1}\right)}^{*}$. So we have the diagonal matrix $\Lambda =Q\left(C_1^{-1}{C}_2^2{\left(C_1^{-1}\right)}^{*}\right)Q^{*}$ where the diagonal elements are eigenvalues of $C_1^{-1}{C}_2^2{\left(C_1^{-1}\right)}^{*}$. Then we obtain, $$\varphi (C_1^2,C_2^2)=\varphi (I,\Lambda )$$

References

1. Perez-Cruz, F. Kullback–Leibler divergence estimation of continuous distributions. Proceedings of the IEEE International Symposium on Information Theory; Toronto, ON, Canada, 6–11 July 2008; pp. 1666-1670.

2. Lee, Y.K.; Park, B.U. Estimation of Kullback–Leibler divergence by local likelihood. Ann. Inst. Stat. Math.; 2006; 58, pp. 327-340. [DOI: https://dx.doi.org/10.1007/s10463-005-0014-8]

3. Anderson, N.; Hall, P.; Titterington, D. Two-sample test statistics for measuring discrepancies between two multivariate probability density functions using kernel-based density estimates. J. Multivar. Anal.; 1994; 50, pp. 41-54. [DOI: https://dx.doi.org/10.1006/jmva.1994.1033]

4. Nguyen, X.; Wainwright, M.J.; Jordan, M.I. Nonparametric estimation of the likelihood ratio and divergence functionals. Proceedings of the IEEE International Symposium on Information Theory; Nice, France, 24–29 June 2007.

5. Wang, Q.; Kulkarni, S.; Verdú, S. Divergence estimation of continuous distributions based on data-dependent partitions. IEEE Trans. Inf. Theory; 2005; 51, pp. 3064-3074. [DOI: https://dx.doi.org/10.1109/TIT.2005.853314]

6. Wang, Q.; Kulkarni, S.; Verdú, S. A nearest-neighbor approach to estimating divergence between continuous random vectors. Proceedings of the IEEE International Symposium on Information Theory; Seattle, WA, USA, 9–14 July 2006; pp. 242-246.

7. Yang, G.L.; Le Cam, L. Asymptotics in Statistics: Some Basic Concepts; Springer: Berlin, Germany, 2000.

8. Kulllback, S.; Leibler, R.A. On information and sufficiency. Ann. Math. Stat.; 1951; 22, pp. 79-86. [DOI: https://dx.doi.org/10.1214/aoms/1177729694]

9. Cichocki, A.; Amari, S. Families of alpha- beta- and gamma- divergences: Flexible and robust measures of similarities. Entropy; 2010; 12, pp. 1532-1568. [DOI: https://dx.doi.org/10.3390/e12061532]

10. Binette, O. A note on reverse Pinsker inequalities. IEEE Trans. Inf. Theory; 2019; 65, pp. 4094-4096. [DOI: https://dx.doi.org/10.1109/TIT.2019.2896192]

11. Sason, I.; Verdú, S. f-divergence inequalities. IEEE Trans. Inf. Theory; 2016; 62, pp. 5973-6006. [DOI: https://dx.doi.org/10.1109/TIT.2016.2603151]

12. Gilardoni, G.L. On Pinsker’s type inequalities and Csiszar’s f-divergences, Part I: Second and Fourth-Order Inequalities. arXiv; 2006; arXiv: cs/0603097

13. Poczos, B.; Schneider, J. On the Estimation of α-Divergences. Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics; Fort Lauderdale, FL, USA, 11–13 April 2011; Volume 15, pp. 609-617.

14. Amari, S. α-Divergence is unique, belonging to both f-divergence and Bregman divergence classes. IEEE Trans. Inf. Theory; 2009; 55, pp. 4925-4931. [DOI: https://dx.doi.org/10.1109/TIT.2009.2030485]

15. Gilberg, A.; Jain, L. If it ain’t broke, don’t fix it: Sparse metric repair. Proceedings of the 55th Annual Allerton Conference on Communication, Control, and Computing (Allerton); Monticello, IL, USA, 3–6 October 2017; pp. 612-619.

16. Ding, R. Visualizing Structures in Financial Time-Series Datasets through Affinity-Based Diffusion Transition Embedding. J. Financ. Data Sci.; 2023; 5, pp. 111-131. [DOI: https://dx.doi.org/10.3905/jfds.2022.1.111]

17. Leonenko, N.N.; Pronzato, L.; Savani, V. A class of Rényi information estimators for multidimensional densities. Ann. Stat.; 2008; 36, pp. 2153-2182. [DOI: https://dx.doi.org/10.1214/07-AOS539]

18. Lehmann, E.L.; Romano, J.P. Testing Statistical Hypotheses; Springer: New York, NY, USA, 2005.