The Ziggurat Algorithm for Random Gaussian Sampling | Hacker News Reader

The Ziggurat Algorithm for Random Gaussian Sampling

heliosphan.org

45 pointsredcalx14 comments

Working through this and other sampling algorithms was a real eye-opening experience for me, and really helped me develop an intuition for statistics in general. It was definitely worthwhile. I'm still rather proud of some code I wrote back in 2009 to implement the Ziggurat algorithm in an unusually general way[1]. It's in Haskell so I apologize in advance for the limited audience, but it's rather well documented (to my surprise... I don't always find that to be the case in my old code, or even my not-so-old stuff).

I always especially liked the fact that it constructed its table lazily and could do so given nothing but a handful of functions and a constant (the PDF itself, its inverse, its definite integral from 0, and the limit of that integral at infinity). From those pieces, it builds the tables and an algorithm to sample, and optionally also recursively transforms those inputs as necessary to lazily construct another ziggurat representation of the tail, ad infinitum. Basically, you can throw any monotone PDF you want at it and not have to worry about solving or approximating the tail analytically.

[1]https://github.com/mokus0/random-fu/blob/master/random-fu/sr...

Very cool. Note that my C# version is fixed to use 128 segments based on George Marsaglia's original calculations; Calculating them dynamically for a given PDF is an impressive achievement indeed. E.g. I did wonder if there was a more optimal number of segments, and of course yours isn't limited to the Gaussian.

The Box–Muller transform method (also mentioned int he article) is kind of fun. You generate two Gaussians at once. You rejection sample from the unit square into the largest contained circle (generates a point about 3/4 of the time). You then roughly log-transform the distance along the radius (sending (0,0) to infinity) and the two points are now independent Gaussians.

Here's my source code for Box-Muller

http://svn.code.sf.net/p/sharpneat/code/branches/V2/src/Shar...

Note it's a lot simpler (more elegant) than the Ziggurat method!

Also note the while loop rejects sqr==0 because that would cause Log(0) to throw an exception. I convinced myself that was correct when I wrote it in 2011, but I'm struggling to remember how I came to that conclusion.

Unfortunately the elegance comes at a cost: Box-Muller is a lot slower due to the sqrt and log. Ziggurat's fast path (which can be hit arbitrarily rarely, by controlling the table size) is just some indexing and comparisons, making it amazingly fast.

Also, Ziggurat can be applied to any distribution with a decreasing PMF (e.g. it's also amazingly fast for the Exponential distribution), whereas Box-Muller relies on special properties of the Normal.

One could say that the high level source code is elegant, but overall its use of log and sqrt invoke a lot of complex machinations (inside the CPU).

Also, IIRC Box-Muller gives incorrect results if you're interested in sampling the tails of the Gaussian.

If you're not interested in the tails, there's also always

(random()+random()+random()+random()-2.0)*sqrt(3.)

as a cheap (in terms of brain power) Gaussian (sigma=1, mean=0) rough approximation :-)

I once used symbolic regression (Eureqa) to generate an approximation to the cumulative distribution. From there you can just draw samples from a uniform distribution, and map them onto a Gaussian distribution. Some software uses an approximated CDF, using polynomials I think.

Here’s a really neat new approach, Olver & Townsend, “Fast inverse transform sampling in one and two dimensions” http://arxiv.org/pdf/1307.1223v1.pdf

Built on top of Chebfun, http://www.chebfun.org

How much faster was it, and how good was the approximation?

The approximation was pretty decent, and it was only a few math functions. I'll try to redo it and compare it to other methods.

Regarding special note 1: If you change the sizes of the segments so that they instead cover equal amounts of probability mass, then you'd be able to sample again from the same segment instead of starting over from the beginning.

But I guess this slows down the fast path a little, so with sufficiently low rejection probability it is not worth it.

Yeh, I think there are probably a few different ways you could 'slice' up the body of the pdff. E.g. I mention using sloped ends to the segments to reduce the amount of pdf area outside of the segments, which would reduce the expensive path by a lot, but it also introduces a division into the fast path. Then again, that division is by a fixed value, so really it's a multiplication (by 1/x), so it might actually be a pretty good solution.

It would be interesting to explore these options further.

Reply on news.ycombinator.com