From marius.hofert at uwaterloo.ca  Wed Mar 25 00:15:11 2015
From: marius.hofert at uwaterloo.ca (Marius Hofert)
Date: Tue, 24 Mar 2015 19:15:11 -0400
Subject: [Rd] Faster version of rgeom()
Message-ID: <CAM3-KjaYnSyNLfdLV6OwfTu-bxEeN_Td3_6od5YacQZuPAzp8w@mail.gmail.com>

Hi,

I recently took a closer look at the implementations of sampling
algorithms for all basic distributions in R. Two seemed inefficient
(slow) to me:
1) rgeom()
2) rexp()
Obviously, Geo(p) and Exp(lambda) have very simple stochastic
representations (floor(Y) for Y ~ Exp(-log(1-p)) and (-log U)/lambda
for U ~ U[0,1]). I thought I could easily beat the (more complicated)
algorithms that R uses for these two distributions. Even Ahrens and
Dieter (1972) (the paper describing the algorithm rexp() uses; see
below) report that the inversion method for sampling Exp(lambda) is
not beatable).

With a student, I put together a small package (attached) comparing
some run times in
the form of a demo (original version R uses vs the above stochastic
representations in both pure R implementations and pure C
implementations). Here are the results (simple screenshots from the
demo attached):
1) rexp: quite different on Ubuntu Linux / OS X. On OS X it seems easy
to beat rexp(), but on Ubuntu not. Still quite a surprise to me... Is
it that log() is implemented differently on different platforms which
was maybe (?) the original reason for using Ahrens and Dieter (1972)
for sampling from Exp(lambda)?
2) rgeom: both platforms more or less show the same results: rgeom()
can be outperformed (by a factor of two or so).

I know that's quite 'philosophical', but I'm wondering whether a
simpler and faster algorithm such as the one-liner for sampling from
Geo(p) is something one could consider to make available in R and keep
rgeom()'s current default for reproducibility (similar to the change
not too long ago where rnorm() was changed to inversion [if I remember
correctly])?

Any comments are welcome. I'm trying to understand and learn with this
(so not complaining that rgeom() is too slow :-) )

Many cheers,
Marius

PS: Ahrens and Dieter (1972): "Computer methods for sampling from the
exponential and normal distributions"




-- 
Marius Hofert, Dr. rer. nat.
Assistant Professor
Department of Statistics and Actuarial Science
Faculty of Mathematics
University of Waterloo
200 University Avenue West, Waterloo, ON, N2L 3G1
+1-519-888-4567, ext. 31394 (office M3 4207)
http://math.uwaterloo.ca/~mhofert
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1_rgeom_comparison_os_x.png
Type: image/png
Size: 397097 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0008.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 1_rgeom_comparison_ubuntu.png
Type: image/png
Size: 63615 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0009.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2_rexp_comparison_os_x.png
Type: image/png
Size: 374402 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0010.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 2_rexp_comparison_ubuntu.png
Type: image/png
Size: 57036 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0011.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3_rgeom_effect_of_prob_os_x.png
Type: image/png
Size: 208246 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0012.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3_rgeom_effect_of_prob_ubuntu.png
Type: image/png
Size: 52621 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0013.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 4_rexp_effect_of_lambda_os_x.png
Type: image/png
Size: 191788 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0014.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 4_rexp_effect_of_lambda_ubuntu.png
Type: image/png
Size: 55383 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0015.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rgeomexp_0.0-1.tar.gz
Type: application/x-gzip
Size: 4388 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20150324/54a9624f/attachment-0001.gz>


[Rd] Faster version of rgeom()
Marius Hofert