Discretization of Kullback-Leibler Divergence

Question

I have a question regarding the Kullback-Leibler divergence. I am working on a project, in which we compare numerically distributions: the empirical distribution of the data compared to an approximate distribution. To measure the divergence, I use the Kullback-Leibler divergence. I've noted that I get negative KL-divergences, which should not be possible (proof here) - so I started looking.

I evaluate both distributions numerically on a grid. So far, I thought that this grid can be arbitrarily designed but I ensure that the maximum and the minimum value of the observed realizations is covered. Since I am also working with heavy-tailed distributions, I wanted to have a grid with more grid elements in the center of the distribution. But as soon as I start to deviate from a linear spaced grid (i.e., meaning the distance between the grid elements is linear), the KL turns negative. I'll append a code in R to replicate the issue. See here:

set.seed(100)

# dimensions
ngrid     = 100  # number of grid elements
ndraws    = 1000 # number of draws
type_grid = "nonlinear"

# create draws from RV
X = rnorm(ndraws, 3, 0.3)
Y = rnorm(ndraws, 3, 0.3)

if(type_grid == "linear"){
  # define grid linear
  xgrid = seq(min(X,Y), max(X,Y), length.out=ngrid) # linear
}else if(type_grid == "nonlinear"){
  # define grid with more mass in the middle
  tmp_dens = density(c(X,Y))
  xgrid = sample(c(X,Y), ngrid, replace=TRUE) + tmp_dens$bw*rnorm(ngrid)
}

# kernel density estimator
fit.dens.X = density(X)
fit.dens.Y = density(Y)

# approximation at grid points
p.X = approx(fit.dens.X$x, fit.dens.X$y, xout = xgrid)$y
p.Y = approx(fit.dens.Y$x, fit.dens.Y$y, xout = xgrid)$y
p.X[is.na(p.X)] = 1e-100
p.Y[is.na(p.Y)] = 1e-100

# kullback-leibler divergence
sum(p.X * (log(p.X/p.Y)))

Any ideas why this is happing or how I can circumvent this? Specifically for heavy-tailed distributions it would be nice to deviate from a linear-spaced grid. Thanks in advance.

Answer 1

In your case, where both distributions $P$ and $Q$ have Lebesgue densities $p$ and $q$, the Kullback-Leibler divergence is given by the integral $$ D_\text{KL}(P \mid Q) = \int_{-\infty}^{\infty} p(x) \log \left( \frac{\ p(x)\ } {q(x)} \right)\, \mathrm{d} x $$ To approximate this integral numerically you cannot simply sum up the values of the integrand, as you seem to do in your code. You have to use quadrature methods from numerical integration, e.g. the Trapezoidal rule. For a non-uniform grid with points $x_0, \ldots, x_N$ this rule gives the formula $$ \int_{a}^{b} f(x)\, \mathrm{d} x \approx \sum_{k=1}^N \frac{f(x_{k-1}) + f(x_k)}{2} \Delta x_k , $$ where $\Delta x_k = x_{k} - x_{k-1}$. In practice you need to choose $a$ and $b$ sufficiently small/large such that the part of the integral outside of $[a,b]$ is negligible.