Evaluating two sets of random samples

Let $p$ be a probability distribution that can be computed tractably for any given point. I use two MCMC methods to generate samples from the distributions. For each MCMC method, I run 1000 Markov chains with random initialization till the chain converges. This gives me two sets of random samples with thousand samples each, one for each MCMC method.

Now I wish to evaluate which set of random samples are more representative of the distribution $p$. In two dimensions, I can do a scatter plot of the points and compare it with a surface plot of the distribution $p$. However, I am not sure what to do in higher dimensions.

Evaluating $\sum_{i=1}^n \log p(x_i)$ for each set individually isn't useful, because a Markov chain that always generates the mode will win. I can compute the derivative of $\sum_{i=1}^n \log p(x_i)$ with respect to the parameters of $p$. The set with the lower norm of derivative is probably a better fit to $p$. However, I am quite sure that there is a much simpler idea. ANy suggestions?

• The Gelman-Rubin statistic is intended for comparing MCMC samples and for providing a rough evaluation as to whether or not they share the same distribution. May 3 '16 at 9:16

You basically want to find the distance between the two distributions, and pick the sampler that has the smaller distance from the probability distribution $p$. The most natural thing to consider is the Total Variation distance since convergence of Markov chains is generally assessed using TV distance. Although, to use this, you might want to make sure you run all the 1000 samples for the same number of iterations before picking the last sample.
The total variation distance between two distributions $P$ and $Q_1$ defined on a space $\mathcal{X}$ is $$\|Q_1 - P\|_{TV} = \sup_{A \subset \mathcal{X}}|Q_1(A) - P(A)|.$$
Intuitively, TV distance is the largest difference between the two distributions over all possible subsets. You can find the TV distance from $P$ for $Q_1$ and $Q_2$, and compare. I believe package distrEx in R can do this for you.
You can also use something like the Kullback-Leibler Divergence. This is more popular in information theory. Note that KL-Divergence is not a distance (which means it is not symmetric). So $KL(P\|Q_1)$ is not the same as $KL(Q_1\|P)$. $KL(P\|Q_1)$ finds the divergence of $Q_1$ from $P$. So you can find $KL(P\|Q_1)$ and $KL(P\|Q_2)$ and compare. In R you can use the package FNN.
• I do not see how to apply this principle in an MCMC setting as the target $P$ is only partly unknown. Do you know of efficient ways of computing the Kullback-Leibler distance in such a case? May 3 '16 at 9:14
• @Xi'an $KL(P\|Q_1) = E_{\pi} [\log(\pi(x))] - E_{\pi}[\log(q(x))]$. So if you were to compare two samplers with the same stationary distribution using KL divergence you could find the log of the empirical density and see which one has the largest $E_{\pi}[\log(q(x)]$. I believe this would require two MCMC runs, one for estimation of $q$ and one for Monte Carlo estimates. May 3 '16 at 16:45