3
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Xi'an
    Commented May 3, 2016 at 9:16

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – Xi'an
    Commented May 3, 2016 at 9:14
  • $\begingroup$ @Xi'an I don't know how I would do this in MCMC. $\endgroup$
    – Greenparker
    Commented May 3, 2016 at 12:11
  • $\begingroup$ @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. $\endgroup$
    – Greenparker
    Commented May 3, 2016 at 16:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.