I'm familiar with the definition of the KL divergence between two discrete distributions
$D_{KL} = D_{KL}\big({\it P}(A) || {\it Q}(B)\big)=\sum_{j=1}^{n} {\it P}(A=a_{j}) \log \Big( \cfrac{{\it P}(A=a_{j})}{{\it Q}(B=b_{j})} \Big)$
and how it approaches the continuous limit as the size of a particular bin width approaches 0.
Implementing this practically with samples from Markov Chain Monte Carlo, I would anticipate that the KL divergence would asymptotically approach the true value as you divide samples into more and more bins.
However, what I find in a simple problem is that the KL divergence continues to drop as the number of bins increase. Below are some figures.
For a parameter $L$, I am comparing an exact distribution vs one I obtained using a Gaussian Process emulator.
Calculating the KL divergence, I find the following behavior when increasing the number of bins:
In the second figure, there are ~50 steps in bin size. My expectation was that the KL divergence would approach a value and then begin to behave erratically as the bin widths became so small that in a given bin, one distribution would have no counts. Why is the KL divergence dropping in such a way and what is the best way to rigorously choose a number of bins for comparing distributions? As far as I've been able to tell, this question hasn't been answered clearly in the literature.
Edit to address a comment from Shimao: If I take samples from two distributions N(0,1) and N(2,1), I get similar behavior.
Edit2: Adding the function I use to calculate the KL Divergence
def KL_divergence(h1, h0, dx):
"""
h0, h1 are histograms
dx is the bin width, which we assume is uniform
Return the information gain of pdf h1 w.r.t. h0, approximated by
riemann sum over histograms of each pdf. Assumes that distributions have the exact same bins
and that bin widths are uniform!
"""
#get nonzero values of posterior, the logarithm is ill-defined for bins with zero
nonzero_idx = [x!=0 for x in h1] # get indices for nonzero bins
h1 = h1[ nonzero_idx ] # filter
h0 = h0[ nonzero_idx ]
nonzero_idx = [x!=0 for x in h0] # get indices for nonzero bins
h1 = h1[ nonzero_idx ] # filter
h0 = h0[ nonzero_idx ]
KLdiv = np.sum( h1 * np.log2(h1/h0) ) * np.prod(dx) # calculate the KL divergence
return KLdiv