On Thu, 24 Dec 2009, Sven Schreiber wrote:
Allin Cottrell schrieb:
> There is a possible catch [in the highspeed ziggurat code]:
> Doornik has an article (google: doornik ziggurat) arguing that
> the standard ziggurat method fails a test for high-dimensional
> randomness (and that one might expect it to fail), and he has
> a workaround which necessarily sacrifices some of the
> efficiency of ziggurat. (Unfortunately his code is free only
> for non-commercial use so it's incompatible with the GPL and
> we can't use it in gretl.)
>
> So far as I can tell Doornik's article hasn't been published and I
> can't find any response from Marsaglia, but I'm taking this
> seriously. Right now I'm running the L'Ecuyer "crush" test on
the
> new gretl code -- should be finished within an hour or two!
It turns out that for the same reason the NumPy people also have not
used the ziggurat so far. Whereas Octave does use it according to its
documentation. Don't know whether they take into account Doornik's
caveat, but it doesn't sound like it and this seems to explain the speed
advantage over Ox (just like in the above new gretl numbers).
What about introducing the (standard, non-Doornik) ziggurat as an
option, configurable through 'set'?
Well, Doornik's specific claim was that the standard
Marsaglia/Tsang ziggurat fails the "collisions" test in the
"crush" test suite. I've now run "crush" on the ziggurat code in
gretl (taken from Jochen Voss, but using the Mersenne Twister as
implemented in glib as the uniform input, rather than gsl's
uniform RNG). It passes all the tests, including "collisions".
We could add a "set" option, but first I think we'd better run
both "crush" and "big crush" (many hours' worth) on both the
ziggurat code and Box-Muller. Right now we don't really have
any basis for thinking that Box-Muller is "safer".
Allin.