12
$\begingroup$

My question is inspired by R's built-in exponential random number generator, the function rexp(). When trying to generate exponentially distributed random numbers, many textbooks recommend the inverse transform method as outlined in this Wikipedia page. I am aware that there are other methods to accomplish this task. In particular, R's source code uses the algorithm outlined in a paper by Ahrens & Dieter (1972).

I have convinced myself that the Ahrens-Dieter (AD) method is correct. Still, I do not see the benefit of using their method compared to the inverse transform (IT) method. AD is not only more complex to implement than IT. There does not seem to be a speed benefit either. Here is my R code to benchmark both methods followed by the results.

invTrans <- function(n)
    -log(runif(n))
print("For the inverse transform:")
print(system.time(invTrans(1e8)))
print("For the Ahrens-Dieter algorithm:")
print(system.time(rexp(1e8)))

Results: 

[1] "For the inverse transform:" 
user     system     elapsed
4.227    0.266      4.597 
[1] "For the Ahrens-Dieter algorithm:"
user     system     elapsed
4.919    0.265      5.213

Comparing the code for the two methods, AD draws at least two uniform random numbers (with the C function unif_rand()) to get one exponential random number. IT only needs one uniform random number. Presumably the R core team decided against implementing IT because it assumed that taking the logarithm may be slower than generating more uniform random numbers. I understand that the speed of taking logarithms can be machine-dependent, but at least for me the opposite is true. Perhaps there are issues around IT’s numerical precision having to do with the singularity of the logarithm at 0? But then, the R source code sexp.c reveals that the implementation of AD also loses some numerical precision because the following portion of C code removes the leading bits from the uniform random number u.

double u = unif_rand();
while(u <= 0. || u >= 1.) u = unif_rand();
for (;;) {
    u += u;
    if (u > 1.)
        break;
    a += q[0];
}
u -= 1.;

u is later recycled as a uniform random number in the remainder of sexp.c. So far, it appears as if

  • IT is easier to code,
  • IT is faster, and
  • both IT and AD possibly lose numerical accuracy.

I would really appreciate if someone could explain why R still implements AD as the only available option for rexp().

Improve this question
$\endgroup$
7
  • 4
    $\begingroup$ With random number generators, "easier to code" is not really a consideration unless you're the one doing it! Speed and accuracy are the only two considerations. (For uniform generators, there's also the period of the generator.) In the old days, A-D was faster. On my Linux box, A-D runs in about 1/2 the time your invTrans function does, and on my laptop in about 2/3 the time. You might want to use microbenchmark for more comprehensive timings, too. $\endgroup$
    – jbowman
    Commented Oct 3, 2016 at 19:39
  • 5
    $\begingroup$ I would suggest we do not migrate it. This seems on-topic to me. $\endgroup$
    – amoeba
    Commented Oct 3, 2016 at 19:40
  • 1
    $\begingroup$ Given that I am unable to think up a single scenario in which rexp(n) would be the bottleneck, the difference in speed is not a strong argument for change (at least to me). I might be more concerned about numeric accuracy, although it's not clear to me which one would be more numerically reliable. $\endgroup$
    – Cliff AB
    Commented Oct 3, 2016 at 19:42
  • 2
    $\begingroup$ @amoeba I think that "What would be the advantages of ..." would be a rephrasing that would be clearly on-topic here, and would not affect any existing answers. I suppose "Why did the people who made R decide to do ..." is really (a) a software-specific question, (b) requires either evidence in the documentation or telepathy, so could arguably be off-topic here. Personally I'd rather the question were rephrased to make it more clearly within the scope of the site, but I don't see this as a strong enough reason to close it. $\endgroup$
    – Silverfish
    Commented Oct 3, 2016 at 20:24
  • 1
    $\begingroup$ @amoeba I've had a go. Not convinced my suggested new title is especially grammatical, and perhaps a few other parts of the question text could do with changing. But I hope this is more clearly on-topic, at least, and I don't think it invalidates or requires changes to either answer. $\endgroup$
    – Silverfish
    Commented Oct 3, 2016 at 20:39

3 Answers 3

Reset to default
9
$\begingroup$

On my computer (pardon my French!):

> print(system.time(rexp(1e8)))
utilisateur     système      écoulé 
      4.617       0.320       4.935 
> print(system.time(rexp(1e8)))
utilisateur     système      écoulé 
      4.589       2.045       6.629 
> print(system.time(-log(runif(1e8))))
utilisateur     système      écoulé 
      7.455       1.080       8.528 
> print(system.time(-log(runif(1e8))))
utilisateur     système      écoulé 
      9.140       1.489      10.623

the inverse transform does worse. But you should watch out for variability. Introducing a rate parameter leads to even more variability for the inverse transform:

> print(system.time(rexp(1e8,rate=.01)))
utilisateur     système      écoulé 
      4.594       0.456       5.047 
> print(system.time(rexp(1e8,rate=.01)))
utilisateur     système      écoulé 
      4.661       1.319       5.976 
> print(system.time(-log(runif(1e8))/.01))
utilisateur     système      écoulé 
     15.675       2.139      17.803 
> print(system.time(-log(runif(1e8))/.01))
utilisateur     système      écoulé 
      7.863       1.122       8.977 
> print(system.time(rexp(1e8,rate=101.01)))
utilisateur     système      écoulé 
      4.610       0.220       4.826 
> print(system.time(rexp(1e8,rate=101.01)))
utilisateur     système      écoulé 
      4.621       0.156       4.774 
> print(system.time(-log(runif(1e8))/101.01))
utilisateur     système      écoulé 
      7.858       0.965       8.819 > 
> print(system.time(-log(runif(1e8))/101.01))
utilisateur     système      écoulé 
     13.924       1.345      15.262

Here are the comparisons using rbenchmark:

> benchmark(x=rexp(1e6,rate=101.01))
  elapsed user.self sys.self
  4.617     4.564    0.056
> benchmark(x=-log(runif(1e6))/101.01)
  elapsed user.self sys.self
  14.749   14.571    0.184
> benchmark(x=rgamma(1e6,shape=1,rate=101.01))
  elapsed user.self sys.self
  14.421   14.362    0.063
> benchmark(x=rexp(1e6,rate=.01))
  elapsed user.self sys.self
  9.414     9.281    0.136
> benchmark(x=-log(runif(1e6))/.01)
  elapsed user.self sys.self
  7.953     7.866    0.092
> benchmark(x=rgamma(1e6,shape=1,rate=.01))
  elapsed user.self sys.self
  26.69    26.649    0.056

So mileage still varies, depending on the scale!

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ On my laptop, the times match with the OP's so closely that I suspect we have the same machine (or at least the same processor). But I think your point here is the speed advantage observed is platform dependent, and given the minimal difference, there's no clear advantage to one over the other in regards to speed. $\endgroup$
    – Cliff AB
    Commented Oct 3, 2016 at 19:52
  • 4
    $\begingroup$ Could you perhaps perform a microbenchmark instead? $\endgroup$
    – Firebug
    Commented Oct 3, 2016 at 20:37
  • 2
    $\begingroup$ The system times appear to measure highly variable overhead, perhaps due to interrupts and memory paging. It is interesting, as @Cliff notes, to see sizable relative differences in performance between systems. On a Xeon with lots of RAM, for instance, I incur almost no system time (0.05 to 0.32 seconds), about 12% longer user time for rexp, 3% shorter user time for -log(runif()), and no variability with a rate parameter ($5.27\pm0.02$ total seconds). We are all implicitly assuming R is achieving times for log and runif comparable to what one would get with a Fortran subroutine. $\endgroup$
    – whuber ♦
    Commented Oct 3, 2016 at 21:10
7
$\begingroup$

This is just quoting the article under the "Algorithm LG: (Logarithm method)" section :

In FORTRAN the algorithm is best programmed directly as $X=-ALOG(REGOL(IR))$ avoiding any subprogram call. The performance was 504 $\mu$sec per sample. Of this time 361 $\mu$sec was taken up by the manufacturer's logarithm routine and 105 $\mu$sec by the generator for REGOL of the uniformly distributed variable $u$. Therefore there was no point in trying to improve speed by writing an assembler function which would have to use the same two subprograms.

So it looks like the authors opted for other methods to avoid this "manufacturer" limitation of slow logarithms. Perhaps this question is then best moved to stackoverflow where someone with knowledge on the guts of R can comment.

Improve this answer
$\endgroup$
6
$\begingroup$

Just running this with microbenchmark; on my machine, R's native approach is uniformly faster:

library(microbenchmark)
microbenchmark(times = 10L,
               R_native = rexp(1e8),
               dir_inv = -log(runif(1e8)))
# Unit: seconds
#      expr      min       lq     mean   median       uq      max neval
#  R_native 3.643980 3.655015 3.687062 3.677351 3.699971 3.783529    10
#   dir_inv 5.780103 5.783707 5.888088 5.912384 5.946964 6.050098    10

For novelty's sake, here's making sure it's not entirely due to having $\lambda = 1$:

lambdas = seq(0, 10, length.out = 25L)[-1L]
png("~/Desktop/micro.png")
matplot(lambdas, 
        ts <- 
          t(sapply(lambdas, function(ll)
            print(microbenchmark(times = 50L,
                                 R_native = rexp(5e5, rate = ll),
                                 dir_inv = -log(runif(5e5))/ll),
                  unit = "relative")[ , "median"])),
        type = "l", lwd = 3L, xlab = expression(lambda),
        ylab = "Relative Timing", lty = 1L,
        col = c("black", "red"), las = 1L,
        main = paste0("Direct Computation of Exponential Variates\n",
                      "vs. R Native Generator (Ahrens-Dieter)"))
text(lambdas[1L], ts[1L, ], c("A-D", "Direct"), pos = 3L)
dev.off()

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.