Kernel density estimation with a Markov chain Monte Carlo sample

Bayesian inference relies on the posterior distribution, which is often estimated with a Markov chain Monte Carlo sampler. The sampler produces a dependent stream of variates from the limiting distribution of the Markov chain, the posterior distribution. When one wishes to display the estimated posterior density, a natural choice is the histogram. However, abundant literature has shown that the kernel density estimator is more accurate than the histogram in terms of mean integrated squared error for an i.i.d. sample. With this as motivation, a kernel density estimation method is proposed that is appropriate for the dependence in the Markov chain Monte Carlo output. To account for the dependence, the cross-validation criterion is modified to select the bandwidth in standard kernel density estimation approaches. A data-driven adjustment to the biased cross-validation method is suggested with introducing the integrated autocorrelation time of the kernel. The convergence of the modified bandwidth to the optimal bandwidth is shown by adapting theorems from the time series literature. Simulation studies show that the proposed method finds the bandwidth close to the optimal value, while standard methods lead to smaller bandwidths under Markov chain samples and hence to undersmoothed density estimates. A study with real data shows that the proposed method has a considerably smaller integrated mean squared error than standard methods. The R package KDEmcmc to implement the suggested algorithm is available on the Comprehensive R Archive Network.