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

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

• Does the Gibbs sampler converge to a global maximum in the presence of multiple modes? For example in case of a Gaussian mixture distribution? – Dhiraj Mar 22 '12 at 19:21

## 2 Answers

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.

• Yup thats what I meant, proportional to, not equals. Anyways, if I'm trying to sample from $P(\theta | D)$, considering that I can calculate something proportional to it, how can I do it that is totally 100% general and not dependent on the particular distributions of $P(\theta)$ and $P(D|\theta)$? – William Aug 16 '11 at 17:56
• I guess you'll have to analytically solve for $P(\theta|D)$, and then run a uniform through the inverse of the cdf unless you do rejection sampling. – highBandWidth Aug 16 '11 at 22:43

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

• can you please give specific clickable references? Not everybody knows how to use Google Scholar, you know. – StasK Aug 20 '11 at 17:21
• @StasK is that ok? – yannick Sep 8 '11 at 7:51
• @StasK if you search for any scientific paper on vanilla Google, then it will bring up a link to some of the Google scholar results at the top (IME) – richiemorrisroe Sep 8 '11 at 11:18
• The slice sampler is very simple to implement, but already exists in a single R file: cs.toronto.edu/~radford/slice.software.html. To get a feel for the issues, you probably should try MH as well. There is no single most efficient algorithm, and more complex approaches will almost certainly perform better. But they are often not worth it. Consider JAGS/Bugs if you're just interesting in quickly sampling from models. – Tristan Sep 8 '11 at 17:41
• @yannick, yep, much better. Name and year are not always meaningful, you know. If you Google naively for something like Cheng et. al. (2010), you may get 400 references, and without good enough knowledge of the specific literature, you would likely be totally lost. – StasK Sep 8 '11 at 20:38