I have lots of experience computing KL divergences for straightforward discrete distributions where I have access to complete probability tables, etc. But I'm a little concerned about my current problem. I have a 2-D distribution, a Gaussian mixture, which I have a contour plot of here:
The above thus represents $p(\theta \mid x_1,\ldots,x_n)$ where $\theta = (\theta_1,\theta_2)$ and the $x_i$s were drawn from a Gaussian mixture model. The modes are at $(0,1)$ and $(1,-1)$.
I have two scatter plots that show MCMC sampling to estimate that posterior distribution.
Both of these attempt to estimate the posterior. However, I would like to have a principled way of telling which distribution is "closer" to the contour plot (the "true" distribution), hence I need the KL divergence.
I think the way I would do this is as follows. For both scatter plots, I would discretize the two dimensions, perhaps at the granularity of 0.05 for both dimensions, and compute the empirical estimate of $p(\theta \mid x_1, \ldots, x_n)$ for each of the discretized $\theta$ values. I think I would round each point in the scatter plot to its closest discretization (e.g., $(0.06,0.09)$ would turn into $(0.05, 0.10)$), and that counts as one spot for that element. Then at the end, after counting up the number of dots in each discretized bin, I'd add perhaps a count of 0.5 to each bin to smooth the distribution and then normalize. Does this make sense?
Oh, I would also have to do this same discretization for that distribution that generated the contour plot, I guess, but that's not so bad because I have the exact function so I can just call it with the bin values.
This is in python so any suggested libraries/methods that could assist in this process would be great.
Remark: I don't think this problem depends on some of the details I said earlier (especially about how $p(\theta \mid x_1, \ldots, x_n)$ gets computed) but I figured I would mention it in case.