Discrete KL Divergence with decreasing bin width 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

 A: The issue that was causing this behavior is that the discrete and continuous versions of the KL divergence approach each other in the limit where the $n_{samples} \rightarrow \infty$ and $width_{bin} \rightarrow 0$.
Holding $n_{samples}$ fixed and increasing $n_{bins}$ as a result is only fulfilling one criteria. The issue of how to choose a number of bins given a number of samples has been answered several times and in introductory texts; Rice (1944), for example, recommends $n_{bins} = 2\sqrt[3]{n_{samples}}$. Implementing this naively (so that sampling fluctuations haven't been taken care of), in the limit $n_{samples} \rightarrow \infty$, the discrete and continuous versions approach agreement for the example suggested by @Shimao in the comment above.
N.B. One must also be careful to compare in the correct units (bits vs. nats).

This convergence can also be seen when calculating KL divergence a set number of times to estimate the sampling error in the histogram below where KL Divergence is along the x axis.

