Каковы преимущества экспоненциального генератора случайных чисел, использующего метод Аренса и Дитера (1972), а не обратным преобразованием?

Мой вопрос вдохновлен встроенным экспоненциальным генератором случайных чисел R , функцией rexp(). При попытке генерировать экспоненциально распределенные случайные числа многие учебники рекомендуют метод обратного преобразования, описанный на этой странице Википедии . Я знаю, что есть другие методы для решения этой задачи. В частности, исходный код R использует алгоритм, изложенный в статье Аренса и Дитера (1972) .

Я убедил себя, что метод Аренса-Дитера (AD) является правильным. Тем не менее, я не вижу преимущества использования их метода по сравнению с методом обратного преобразования (IT). AD не только сложнее в реализации, чем IT. Похоже, выигрыша в скорости тоже нет. Вот мой код R для сравнения обоих методов с последующими результатами.

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)))

Полученные результаты:

[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

Сравнивая код для двух методов, AD рисует как минимум два одинаковых случайных числа (с функцией Cunif_rand() ), чтобы получить одно экспоненциальное случайное число. Этому нужно только одно равномерное случайное число. Предположительно, основная команда R отказалась от внедрения IT, поскольку предположила, что логарифм может быть медленнее, чем генерация более однородных случайных чисел. Я понимаю, что скорость взятия логарифмов может зависеть от машины, но, по крайней мере, для меня все наоборот. Возможно, есть проблемы, связанные с числовой точностью ИТ, связанные с сингулярностью логарифма при 0? Но тогда исходный код R sexp.cпоказывает, что реализация AD также теряет некоторую числовую точность, потому что следующая часть кода C удаляет начальные биты из единого случайного числа 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 позже переработаны в виде однородного случайного числа в остальной части sexp.c . Пока что кажется, что

• ЭТО легче кодировать,
• Это быстрее, и
• И IT, и AD, возможно, теряют числовую точность.

Я был бы очень признателен, если бы кто-то мог объяснить, почему R по-прежнему использует AD как единственный доступный вариант rexp().

На моем компьютере (простите мой французский!):

> 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

обратное преобразование делает хуже. Но вы должны следить за изменчивостью. Введение параметра скорости приводит к еще большей изменчивости для обратного преобразования:

> 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 

Вот сравнение с использованием 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

Таким образом, пробег все еще варьируется, в зависимости от масштаба!

Это просто цитирование статьи в разделе «Алгоритм LG: (метод логарифма)»:

$X=-ALOG(REGOL(IR))$$\mu$$\mu$$\mu$$u$

Таким образом, похоже, что авторы выбрали другие методы, чтобы избежать этого «изготовителя» ограничения медленных логарифмов. Возможно, тогда этот вопрос лучше всего перенести в stackoverflow, где кто-то, обладающий знаниями в области R, может прокомментировать.

Просто запускаю это с microbenchmark; на моей машине родной подход R одинаково быстрее:

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

$\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()