5:16 p.m.
Hello all,
Some news and a warning.
For some time now we've been using the Marsaglia-Tsang "Ziggurat"
algorithm for generating (pseudo-)random normal variates. Recently,
Jack noticed that we were considerably slower than Octave at doing
this, so we compared the respective codes.
It turned out that the Octave Ziggurat code (written by David Bateman)
was indeed better than what we had: not only faster, but also gives
denser coverage of the real line. We checked it using l'Ecuyer's "Big
Crush" suite and it passed with flying colors. Therefore in current
git we have switched to Bateman's implementation.
This is a nice improvement, but here's the warning: unfortunately, it
is impossible to preserve the sequence of pseudo-random values
produced for any given seed, across the change in algorithm. So in
current git (and the upcoming 2016d release) you will not get the same
sequences as with gretl 2016c and earlier (back to gretl 1.8.7, when
we last changed the code for generating normals).
If you have Monte Carlo code that runs many replications, results
should be qualitatively very similar, but "single-case" simulations
may produce quite different output for a given seed.
--
Allin Cottrell
Department of Economics
Wake Forest University, NC
8:14 p.m.
Am 15.11.2016 um 18:16 schrieb Allin Cottrell:
Hello all,
Some news and a warning.
For some time now we've been using the Marsaglia-Tsang "Ziggurat"
algorithm for generating (pseudo-)random normal variates. Recently, Jack
noticed that we were considerably slower than Octave at doing this, so
we compared the respective codes.
It turned out that the Octave Ziggurat code (written by David Bateman)
was indeed better than what we had: not only faster, but also gives
denser coverage of the real line. We checked it using l'Ecuyer's "Big
Crush" suite and it passed with flying colors. Therefore in current git
we have switched to Bateman's implementation.
Thanks -- could you perhaps also comment on the comparison in this (not
brand new) talk by Dirk Eddelbüttel:
http://dirk.eddelbuettel.com/papers/useR2014_rcppziggurat.pdf
Check out slide number 11 (or 13, depending on the counting), where
"ZiggGretl" appears.
cheers,
sven
12:03 a.m.
On Tue, 15 Nov 2016, Sven Schreiber wrote:
Am 15.11.2016 um 18:16 schrieb Allin Cottrell:
> Hello all,
>
> Some news and a warning.
>
> For some time now we've been using the Marsaglia-Tsang "Ziggurat"
> algorithm for generating (pseudo-)random normal variates. Recently, Jack
> noticed that we were considerably slower than Octave at doing this, so
> we compared the respective codes.
>
> It turned out that the Octave Ziggurat code (written by David Bateman)
> was indeed better than what we had: not only faster, but also gives
> denser coverage of the real line. We checked it using l'Ecuyer's "Big
> Crush" suite and it passed with flying colors. Therefore in current git
> we have switched to Bateman's implementation.
Thanks -- could you perhaps also comment on the comparison in this (not brand
new) talk by Dirk Eddelbüttel:
http://dirk.eddelbuettel.com/papers/useR2014_rcppziggurat.pdf
Check out slide number 11 (or 13, depending on the counting), where
"ZiggGretl" appears.
I don't have much to say about that reference since I'm not sure
from looking at Dirk's slides how the timings were obtained. That
said, I can say that my timings for 200 replications of (1 million
random normals) show R's rnorm() function -- from the current R
release -- as taking about 5 times as long as gretl and Octave.
Not in the ballpark.
As regards gretl vs Octave, we're now a little faster (using
basically the same normals code, but gretl's generation of the
underlying uniform variates is faster).
Maybe Dirk is using a "souped up" alternative to R's standard
rnorm(). That wasn't clear to me.
Allin
12:41 a.m.
On Tue, 15 Nov 2016, Allin Cottrell wrote:
I don't have much to say about that reference since I'm not sure from looking
at Dirk's slides how the timings were obtained. That said, I can say that my
timings for 200 replications of (1 million random normals) show R's rnorm()
function -- from the current R release -- as taking about 5 times as long as
gretl and Octave. Not in the ballpark.
As regards gretl vs Octave, we're now a little faster (using basically the
same normals code, but gretl's generation of the underlying uniform variates
is faster).
Maybe Dirk is using a "souped up" alternative to R's standard rnorm(). That
wasn't clear to me.
I ran a few normality tests by following Dirk's recipes. Not as many as he
did, but I'd be inclined to say that we are ok, and by a very reassuring
margin too.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:28 a.m.
Am 16.11.2016 um 01:03 schrieb Allin Cottrell:
I don't have much to say about that reference since I'm not sure from
looking at Dirk's slides how the timings were obtained. That said, I can
say that my timings for 200 replications of (1 million random normals)
show R's rnorm() function -- from the current R release -- as taking
about 5 times as long as gretl and Octave. Not in the ballpark.
That would seem to be compatible with the comparison of R's "BM" on a
previous slide with gretl's Ziggurat, roughly 350ms vs. 70ms. (I think
Box-Muller is the default on R, at least in the past.)
Maybe Dirk is using a "souped up" alternative to R's standard rnorm().
That wasn't clear to me.
My impression is that it's an additional package ("RcppZiggurat"?),
which then perhaps provides this function zrnorm() on slide 22.
The difference between "Zigg" and "ZiggGretl" on hte slide 13 is about
20ms vs. 70ms. When you wrote "considerably slower than Octave", what
was the factor you had found?
thanks,
sven
8:31 p.m.
On Wed, 16 Nov 2016, Sven Schreiber wrote:
Am 16.11.2016 um 01:03 schrieb Allin Cottrell:
>
> I don't have much to say about that reference since I'm not sure from
> looking at Dirk's slides how the timings were obtained. That said, I can
> say that my timings for 200 replications of (1 million random normals)
> show R's rnorm() function -- from the current R release -- as taking
> about 5 times as long as gretl and Octave. Not in the ballpark.
That would seem to be compatible with the comparison of R's "BM" on a
previous slide with gretl's Ziggurat, roughly 350ms vs. 70ms. (I think
Box-Muller is the default on R, at least in the past.)
>
> Maybe Dirk is using a "souped up" alternative to R's standard rnorm().
> That wasn't clear to me.
My impression is that it's an additional package ("RcppZiggurat"?), which
then perhaps provides this function zrnorm() on slide 22.
You're right. Sorry, that should have been obvious to me.
The difference between "Zigg" and "ZiggGretl" on
hte slide 13 is about 20ms
vs. 70ms. When you wrote "considerably slower than Octave", what was the
factor you had found?
That seems to depend quite substantially on hardware. When Jack
first raised the normals speed issue, he showed gretl taking more
than twice as long as octave on his laptop. When I tried the same
experiment on my machines, I found a difference of around 20 or 30
percent.
Allin
4:23 p.m.
On Wed, 16 Nov 2016, Sven Schreiber wrote:
My impression is that it's an additional package
("RcppZiggurat"?),
which then perhaps provides this function zrnorm() on slide 22.
I've now looked into this in detail.
The speed with which a program generates normal variates depends on
two factors: the speed of the underlying uniform generator and the
speed of conversion to normal.
In Jack's tests of gretl against Octave -- which set us off on
revising the gretl RNG -- it transpired that gretl's uniform generator
(SFMT) was faster than Octave's but our Ziggurat conversion was
slower, making us somewhat slower overall for normals. Now that we've
switched to Bateman's Ziggurat implementation we're a bit faster
overall.
Now for RcppZiggurat: I've experimented with two of Dirk's functions,
zrnorm() and zrnormR(). The first of these is overall very fast, due
to the fact that the uniform generator he's using, SHR3, is very fast.
However, it has a rather short period by today's standards, 2^32-1 as
against 2^19937-1 for the MT variant used by gretl. The second
function uses the same Ziggurat code but employs R's standard MT
uniform generator, which makes it a bit slower.
Here's what I'm seeing for timings of 200 replications of generating 1
million normals (Nehalem CPU):
gretl: 3.978 seconds
Octave: 5.185
R, zrnorm: 3.275
R, zrnormR: 4.995
So in this comparison class we have the best time for MT + Ziggurat,
though RcppZiggurat is faster with SHR3 + Ziggurat.
Allin
5:15 p.m.
Am 17.11.2016 um 17:23 schrieb Allin Cottrell:
Here's what I'm seeing for timings of 200 replications of
generating 1
million normals (Nehalem CPU):
gretl: 3.978 seconds
Octave: 5.185
R, zrnorm: 3.275
R, zrnormR: 4.995
So in this comparison class we have the best time for MT + Ziggurat,
though RcppZiggurat is faster with SHR3 + Ziggurat.
Thanks for the interesting facts. So we can only speculate why on Dirk's
slides "ZiggGretl" is so much slower. (Perhaps that really only tested
the Ziggurat step of the algorithm, not the overall RNG performance.)
thanks,
sven
2962
days inactive
2964
days old
7 comments
3 participants
participants (3)
-
Allin Cottrell -
Riccardo (Jack) Lucchetti -
Sven Schreiber