Simulating KL Divergence between Cauchy RV and the MLE estimate of the RV - Multimodality seems wrong

Question

I'm working on a (what I think is a fairly simple, straightforward) explanation of how it's really hard to approximate distributions with fat tails accurately in the tails. I started looking at an example, (which ends up harder than I expected,) for KL-divergence between a Cauchy RV and an approximation of the Cauchy RV using MLE. I realized I'm not sure how to do this analytically, and it seems really messy. (Any suggestions for how to do this would be appreciated - but that's a more math-oriented problem. If you want to help me, see question here: https://math.stackexchange.com/questions/2733994/finding-kl-divergence-between-a-distribution-and-the-mle-estimate-of-that-distri

In any case, to try to get a handle on this, however, I built a simulation in R, and the histogram of the KL divergences seems really really strange, and I think there is something wrong, or I'm creating an artifact because of how I'm simulating. (I have tried changing the number of breaks for the histogram, and the number of samples. It doesn't help.)

So first, I wrote a function to simulate a single KL-Divergence

simulate_cauchy_KLD <- function() {
a     <- rcauchy(5000, 0, 1)
coefp <- coef(fitdist(a, 'cauchy', start = list(location = 0, scale = 1)))
h  <- hist(a, breaks = 10000,plot=F) # Only needed if plot=T --- freq = F, ylim = c(0,.35), xlim = c(-10, 10)
ppar<- dcauchy(h$mids, coefp[1], coefp[2])
trus<- h$density
rm(h)
ppar <- ppar[which(trus>0)]
trus<- trus[which(trus>0)]
return(sum(trus*log10(trus/ppar)))

}

Now, I simulate it a bunch of times.

trials=3000 #Takes 1 minute on my system. Seems like a good number.
KLD_cauchy = rep(0,trials)

for(i in 1:trials){
KLD_cauchy[i]<- simulate_cauchy_KLD()
}
hist(KLD_cauchy)

I get a Histogram with these really unexpected groupings:

I'm currently at a loss to understand this. I'd really like either an explanation of what I'm doing wrong, or some idea why this would occur.

After discussion, it was suggested that I try the same process with a normal distribution, and I got the below.

Answer 1

There are several Monte Carlo difficulties with this approach:

1. The first one is that, when handling a sample $x_1,\ldots,x_T$ from a density $f$, the Monte Carlo average $$\frac{1}{T}\sum_{t=1}^T f(x_i) \log \{f(x_i)/g(x_i)\}$$ converges to$$\int f(x) \log \{f(x)\big/g(x)\} f(x)\,\text{d}x=\int \log \{f(x)\big/g(x)\} f^2(x)\,\text{d}x$$ not to $$\int \log \{f(x)\big/g(x)\} f(x)\,\text{d}x$$
2. The second one is that the sample $x_1,\ldots,x_T$ used in the R code is Cauchy and not generated from the histogram distribution $h_T$, which means that the Monte Carlo average $$\frac{1}{T}\sum_{t=1}^T h_T(x_i) \log \{h_T(x_i)/g(x_i)\}$$ converges to$$\int h_T(x) \log \{h_T(x)\big/g(x)\} f(x)\,\text{d}x$$
3. The third one is that even if the above Kullback-Leibler integral should be well-defined (since the histogram distribution $h_T$ has a finite support) it should vary considerably from Cauchy sample to Cauchy sample. For instance, the $1/5000$th quantile of the Cauchy is -1591... This does not explain for the multimodality but certainly for the instability.

A final point is to wonder about the utility of deriving this Kullback-Leibler distance. Is it to compare MLE with other estimators?

Answer 2

The well-separated bumps in your simulation histogram are a consequence of the argument breaks=10000. Given such a request, R will create roughly 10000 equal-width bins, but the bin width, while prettified, may vary from run to run, depending on the presence of extreme outliers in your sample. The presence of outliers has several consequences:

1. To handle an outlier, the requested number of bins must cover a larger region, forcing a larger bin width.

2. You've used density values, i.e., the height at each bin, in your KLD calculation. You should be using the probability of each bin (i.e. height of bin times width of bin), not the height alone. By omitting what amounts to a normalizing factor, the KLD scales differently with different bin widths. If the bin width were the same from run to run, omitting bin width would not be an issue.

3. An extreme outlier contributes a large value of trus/ppar to your KLD. Thus a large bin width correlates not only with the presence of an outlier, but also with with a high KLD.

4. Your line ppar <- ppar[which(trus>0)] has the effect of deleting mass from the "true" density ppar. When bin width is large, more mass gets deleted than with a smaller bin width. This loss of mass appears to have a non-negligible effect on KLD, as you can see by plotting KLD vs bin width or KLD vs sum(ppar) after deletion of mass.

One fix would be to multiply bin height by bin width in the KLD calculation. The widely separated bumps go away when you do this, but the resulting KLD histogram can still exhibit bimodality because of the loss of mass caused by (4).

In my simulations if I force a bin width consistent across all trials (e.g. by using breaks=seq(-10,10,.01)) or if I renormalize ppar, I find the resulting KLD much better behaved. (There is some difficulty in forcing a fixed binwidth when the distribution is heavy-tailed.)

Answer 3

I am still unsure why the simulation is doing this, but at the prompting of @whuber, I tried simulating using R's 'integrate' function using an approach suggested in this question; Estimate the Kullback-Leibler divergence. This removes the multimodality. (I am still unsure exactly what created this artifact.)

integrate_cauchy_KLD <- function() {
  Samples = rcauchy(5000, 0, 1)
  coefp <- coef(fitdist(Samples, 'cauchy', start = list(location = 0, scale = 1)))
  KL_Integrand = function(x){dcauchy(x) * (log2(dcauchy(x))- log2(dcauchy(x, location = coefp[1], scale = coefp[2])))}
  return(integrate(KL_Integrand, -Inf, Inf))

}

Simulating as before, and generating the histogram;

trials=3000
KLD_cauchy_int = rep(0,trials)

for(i in 1:trials){
  KLD_cauchy_int[i]<- integrate_cauchy_KLD()[1]$value
}
hist(KLD_cauchy_int, breaks = 1000)