Kernel density-based Kullback-Leibler divergence estimation in any dimension

Disclaimer: this function doesn't use binning and/or the fast Fourier transform and hence, it is extremely slow even for moderate datasets. For this reason, it is not exported currently.

Usage

kld_est_kde(X, Y, hX = NULL, hY = NULL, rule = c("Silverman", "Scott"))

Arguments

X, Y: n-by-d and m-by-d matrices, representing n samples from the true distribution $P$ and m samples from the approximate distribution $Q$, both in d dimensions. Vector input is treated as a column matrix.
hX, hY: Positive scalars or length d vectors, representing bandwidth parameters (possibly different in each component) for the density estimates of $P$ and $Q$, respectively. If unspecified, a heurestic specified via the rule argument is used.
rule: A heuristic for computing arguments hX and/or hY. The default "silverman" is Silverman's rule $$h_i = \sigma_i\left(\frac{4}{(2+d)n}\right)^{1/(d+4)}.$$ As an alternative, Scott's rule "scott" can be used, $$h_i = \frac{\sigma_i}{n^{1/(d+4)}}.$$

Value

A scalar, the estimated Kullback-Leibler divergence $\hat D_{KL}(P||Q)$.

Details

This estimation method approximates the densities of the unknown distributions $P$ and $Q$ by kernel density estimates, using a sample size- and dimension-dependent bandwidth parameter and a Gaussian kernel. It works for any number of dimensions but is very slow.

Examples

# KL-D between two samples from 1-D Gaussians:
set.seed(0)
X <- rnorm(100)
Y <- rnorm(100, mean = 1, sd = 2)
kld_gaussian(mu1 = 0, sigma1 = 1, mu2 = 1, sigma2 = 2^2)
#> [1] 0.4431472
kld_est_kde1(X, Y)
#> [1] 0.3773503
kld_est_nn(X, Y)
#> [1] 0.2169136
kld_est_brnn(X, Y)
#> [1] 0.1865685

# KL-D between two samples from 2-D Gaussians:
set.seed(0)
X1 <- rnorm(100)
X2 <- rnorm(100)
Y1 <- rnorm(100)
Y2 <- Y1 + rnorm(100)
X <- cbind(X1,X2)
Y <- cbind(Y1,Y2)
kld_gaussian(mu1 = rep(0,2), sigma1 = diag(2),
             mu2 = rep(0,2), sigma2 = matrix(c(1,1,1,2),nrow=2))
#> [1] 0.5
kld_est_kde2(X, Y)
#> [1] 0.07137239
kld_est_nn(X, Y)
#> [1] 0.1841371
kld_est_brnn(X, Y)
#> [1] 0.1854864