Jakie są zalety wykładniczego generatora losowego wykorzystującego metodę Ahrensa i Dietera (1972) zamiast transformacji odwrotnej?

Moje pytanie jest inspirowane wbudowanym generatorem wykładniczej liczby losowej R. , funkcją rexp(). Podczas próby generowania wykładniczych liczb losowych rozkładanych wykładniczo wiele podręczników zaleca metodę transformacji odwrotnej opisaną na tej stronie Wikipedii . Wiem, że istnieją inne metody realizacji tego zadania. W szczególności kod źródłowy R korzysta z algorytmu przedstawionego w artykule Ahrensa i Dietera (1972) .

Przekonałem się, że metoda Ahrensa-Dietera (AD) jest poprawna. Mimo to nie widzę korzyści z zastosowania ich metody w porównaniu z metodą transformacji odwrotnej (IT). AD jest nie tylko bardziej skomplikowany do wdrożenia niż IT. Wydaje się, że nie ma również korzyści z prędkości. Oto mój kod R do porównania obu metod, a następnie wyników.

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

Wyniki:

[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

Porównując kod dla dwóch metod, AD rysuje co najmniej dwie jednolite liczby losowe (z funkcją Cunif_rand() ), aby uzyskać jedną wykładniczą liczbę losową. IT potrzebuje tylko jednej jednolitej liczby losowej. Prawdopodobnie główny zespół R zdecydował się nie wdrażać technologii informatycznych, ponieważ zakładał, że przyjęcie logarytmu może być wolniejsze niż generowanie bardziej jednolitych liczb losowych. Rozumiem, że szybkość przyjmowania logarytmów może zależeć od maszyny, ale przynajmniej dla mnie jest odwrotnie. Być może istnieją problemy związane z dokładnością liczbową IT związaną z osobliwością logarytmu przy 0? Ale potem kod źródłowy R sexp.cujawnia, że ​​implementacja AD również traci pewną dokładność liczbową, ponieważ następna część kodu C usuwa wiodące bity z jednolitej liczby losowej 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 jest następnie zawracana jako jednolity losowej liczby w pozostałej sexp.c . Jak dotąd wydaje się, jakby

• Łatwiej jest kodować,
• IT jest szybsze i
• zarówno IT, jak i AD prawdopodobnie tracą dokładność liczbową.

Byłbym naprawdę wdzięczny, gdyby ktoś mógł wyjaśnić, dlaczego R nadal implementuje AD jako jedyną dostępną opcję dla rexp().

Na moim komputerze (wybacz mój francuski!):

> 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

odwrotna transformacja gorzej. Ale powinieneś uważać na zmienność. Wprowadzenie parametru szybkości prowadzi do jeszcze większej zmienności transformacji odwrotnej:

> 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 

Oto porównania przy użyciu 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

Tak więc przebieg nadal jest różny, w zależności od skali!

Cytuję tylko ten artykuł w sekcji „Algorytm LG: (metoda logarytmiczna)”:

W FORTRAN algorytm najlepiej zaprogramować bezpośrednio jako unikając jakiegokolwiek wywołania podprogramu. Wydajność wynosiła 504 sec na próbkę. Ten czas 361 s rozpuszcza się za pomocą rutynowych logarytm producenta i 105 s przez generator do REGOL z równomiernie zmiennej . Dlatego nie było sensu próbować poprawiać prędkości, pisząc funkcję asemblera, która musiałaby używać tych samych dwóch podprogramów.μ μ μ $X=-ALOG(REGOL(IR))$$\mu$$\mu$$\mu$$u$

Wygląda więc na to, że autorzy wybrali inne metody, aby uniknąć ograniczenia „wolnych producentów” w logarytmach. Być może to pytanie najlepiej przenieść na stackoverflow, w którym ktoś z wiedzą na temat odwagi R może komentować.

Po prostu uruchamiam to z microbenchmark; na moim komputerze natywne podejście R jest jednolicie szybsze:

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

Na wszelki wypadek, upewnijmy się, że nie jest to całkowicie spowodowane tym, że :$\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()