What is the most computationally efficient way to sample from an unnormalized density?

Question

EDIT: After doing some more research, it seems slice sampling could be the way to go. I've seen it mentioned a lot in the context of how to sample the univariate distributions required for Gibbs Sampling. Any thoughts/opinions on this? Any other algorithms that might work better?

This is a half statistical and half methodological question, so please keep an open mind as to any type of solution.

I've been using Gibbs Sampling a lot recently to run Bayesian models and to sample from their posterior predictive distributions. Typically, I've been using the conjugate distributions, e.g. a linear regression with Gaussian-distributed $\beta$ with a Gaussian prior, and a Gamma-distributed noise term $\epsilon_t$. But it occurred to me that part of the real value of MCMC is supposed to be that, unlike message-passing and variational algorithms, the distributions don't necessarily need to be conjugate. Theoretically, I should be able to use any distribution, and the Gibbs Sampling should still work.

So, once I have the distribution $P(\theta | D) = P(D | \theta)P(\theta)$, what I have previously done is work out the form of the distribution on paper, transform it into a distribution I know about and figure out its parameters so I can sample from it.

My questions are:

1. In situations like this is it better to use regular Metropolis-Hastings rather than Gibbs Sampling?
2. Do you know of a good library for doing this kind of symbolic manipulation, so I don't have to do it by hand?
3. If a conjugate distribution can't be found, for example if I'm using some crazy distribution I invented myself, what algorithms can most efficiently sample from $P(\theta | D)$ despite the fact that I don't know its normalization constant?

My main interest is in #3, since it's sort of a catch-all solution -- I'm aware of the family of rejection sampling algorithms, but it seems crazy that while running a sampling algorithm I'd have to spawn yet another sampling algorithm. What is the most efficient way to go about this?

Thanks very much!

Jason

Answer 1

First of all, $P(\theta,D)=P(D|\theta)P(\theta)$ and not $P(\theta|D)$. Perhaps it is a type since you refer to it as an un-normalized version. Secondly, you may not need to run two rejection sampling algorithms since the prior $P(\theta)$ can usually be sampled directly and then you can reject it with $P(D|\theta)$. What is the dimension of $\theta$? Only if you are interested in the joint (or marginals, expectations etc.) of a many dimensional distribution does Gibbs make sense. I think Metropolis Hastings might be beneficial if direct sampling from $P(\theta)$ and rejections using $P(D|\theta)$ leads to a very low overall acceptance rate.

There are some symbolic math capabilities in Mathematica and SymPy that I know of.

Answer 2

It all depends on three aspects of your system:

• the number of modes
• dimensionality
• the correlation structure of $\theta$

If you expect all parameters to have a unique solution, and your posterior to be unimodal, sampling is quite easy with all the methods you just cited, and I wouldn't really bother looking further, since you will always end up in the global maximum of the posterior.

If $\theta$ has lots of coordinates, slice sampling might still work, although I don't know that method really well. However, being a form of Gibbs sampling, I have the feeling that Slice sampling has the following drawback. Reading Neal's 2003 article page 712, slice sampling is performed in the following way

(a) Draw a real value, $y$, uniformly from $(0,f(x0))$, thereby defining a horizontal “slice”: $S = {x : y < f (x)}$. Note that $x_0$ is always within $S$.

(b) Find an interval, $I = (L, R)$, around $x_0$ that contains all, or much, of the slice.

(c) Draw the new point, $x_1$, from the part of the slice within this interval.

If I'm not mistaken, step (b) could actually be very painful if the coordinates of $\theta$ (in Neal's notation, the $x_i$) are highly correlated. You would end up having an infinitely small patch $I$, and convergence would be slow.

I am personnally a big fan of Hybrid monte carlo (Duane 1987), which is a monte carlo combined with molecular dynamics (e.g. overrelaxation, see also the end of Neal's paper). It has the advantage that you propose concerted changes of $\theta$'s coordinates, since they are modified at the same time. It comes with additional parameters you need to tune, but I believe it's really powerful.

EDIT: Here are the references. Since it's copyrighted material, I link to the journal pages where you can try to download the papers. This one is free This one needs a paid subscription, but you can find it on google if you look for it