I've noticed that a lot of physics Monte Carlo simulations make extensive use of rejection sampling, rather than inverse transform sampling. In my research, I'm sampling random energy transfers from energy loss distributions, ie: f(E) represents the relative probability of losing E energy in a collision between two particles. So I need to be able to generate random energy transfers such that the distribution of said transfers will converge to f(E).

I personally use inverse transform sampling all the time, and I'm wondering why people use rejection sampling.

Rejection sampling requires that you generate potentially many random numbers, whereas inverse transform sampling gives you a random variate from the desired distribution every time you generate your uniform variate. Rejection sampling also requires that you can test your variate selection using the original f(E), which could be a complicated function or a large numerical dataset in its own right.

For inverse transform sampling, I generate lookup tables: if I generate a uniform random variate R, I find the closest value in my lookup table and select my non-uniform random variate, which is a O(1) operation, since it's just indexing into an array, so run-time speed is not an issue.

Wikipedia mentions that generating those lookup tables may be computationally difficult, but the example they give of selecting from a random normal distribution is pretty trivial to numerically integrate, so I could generate the probability distribution quite readily. And if something has discontinuities that make it hard to integrate, it may not be any easier to generate a bounding box that's suitable for rejection sampling.

Are there any concrete examples of where this won't work? I suppose for many-dimensional large datasets, the lookup table might become prohibitively large, although I've happily used it with 3D probability distributions with no problems. And if you were generating a new probability distribution on the fly, so that you had to regenerate the lookup tables, then rejection sampling would be the way to go, but that's not the case in most of the papers I've been reading.

  • $\begingroup$ I imagine simplicity. $\endgroup$ – copper.hat Jun 3 '15 at 23:57
  • $\begingroup$ The Metropolis algorithm wiki en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm mentions one of its first usages was for the canonical ensemble which can very difficult constants normalization factors to calculate. $\endgroup$ – muaddib Jun 4 '15 at 0:31
  • $\begingroup$ Is simplicity really a good enough reason in this age of software libraries and packages? And the Metropolis-Hastings algorithm is another good alternative to rejection sampling, although I suppose it could be considered a very complicated form of rejection sampling too. $\endgroup$ – Amanda Ellaway Jun 4 '15 at 2:08

If your function is simple then you can use the "inverse transform sampling". But if your function is complicated - it is very difficult to calculate it's inverse transform of your function then you should consider using the rejection sampling algorithm or other alternatives.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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