

The Ziggurat Algorithm for Random Gaussian Sampling - redcalx
http://heliosphan.org/zigguratalgorithm/zigguratalgorithm.html?hn

======
mokus
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...](https://github.com/mokus0/random-fu/blob/master/random-
fu/src/Data/Random/Distribution/Ziggurat.hs)

~~~
redcalx
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.

------
jmount
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.

~~~
redcalx
Here's my source code for Box-Muller

[http://svn.code.sf.net/p/sharpneat/code/branches/V2/src/Shar...](http://svn.code.sf.net/p/sharpneat/code/branches/V2/src/SharpNeatLib/Utility/BoxMullerGaussianSampler.cs)

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.

~~~
dbaupp
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.

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

~~~
onnoonno
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 :-)

------
Houshalter
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.

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

~~~
Houshalter
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.

------
im3w1l
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.

~~~
redcalx
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.

