Allin Cottrell schrieb:
On Wed, 23 Dec 2009, Sven Schreiber wrote:
> Talha Yalta schrieb:
>>> ...I have "borrowed" the code Jochen
>>> Voss wrote for gsl to implement the ziggurat method for generating
>>> random normal values. Here's what I'm seeing after making this
>>> change:
>>>
>>> loop vectorized
>>> octave 2.59 0.53
>>> gretl 0.87 0.39
>>> ox 0.85 0.56
>> Wow!
>>
> Yes I was gonna say that as well!
There is a possible catch: 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'?
cheers,
sven