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?
 A: 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.
