Was sind die Vorteile eines exponentiellen Zufallsgenerators nach der Methode von Ahrens und Dieter (1972) anstelle einer inversen Transformation?

Meine Frage ist inspiriert R ‚s built-in exponentiellem Zufallszahlengenerator, die Funktion rexp(). Bei dem Versuch, exponentiell verteilte Zufallszahlen zu generieren, empfehlen viele Lehrbücher die inverse Transformationsmethode, wie auf dieser Wikipedia-Seite beschrieben . Mir ist bewusst, dass es andere Methoden gibt, um diese Aufgabe zu erfüllen. Insbesondere verwendet der Quellcode von R den Algorithmus, der in einem Artikel von Ahrens & Dieter (1972) beschrieben wurde .

Ich habe mich davon überzeugt, dass die Ahrens-Dieter (AD) -Methode richtig ist. Dennoch sehe ich keinen Vorteil in der Verwendung ihrer Methode im Vergleich zur inversen Transformationsmethode (IT). AD ist nicht nur komplexer zu implementieren als IT. Es scheint auch keinen Geschwindigkeitsvorteil zu geben. Hier ist mein R- Code zum Benchmarking beider Methoden, gefolgt von den Ergebnissen.

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

Ergebnisse:

[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

Beim Vergleich des Codes für die beiden Methoden zeichnet AD mindestens zwei einheitliche Zufallszahlen (mit der C- Funktion unif_rand()), um eine exponentielle Zufallszahl zu erhalten. Die IT benötigt nur eine einheitliche Zufallszahl. Vermutlich hat sich das R- Kernteam gegen die Implementierung der IT entschieden, da davon ausgegangen wurde, dass der Logarithmus langsamer ist als die Generierung einheitlicherer Zufallszahlen. Ich verstehe, dass die Geschwindigkeit der Logarithmen maschinenabhängig sein kann, aber zumindest für mich ist das Gegenteil der Fall. Vielleicht gibt es Probleme mit der numerischen Genauigkeit der IT, die mit der Singularität des Logarithmus bei 0 zu tun haben? Aber dann der R- Quellcode sexp.czeigt, dass die Implementierung von AD auch eine gewisse numerische Genauigkeit verliert, da der folgende Teil des C-Codes die führenden Bits aus der einheitlichen Zufallszahl u entfernt .

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 wird später im Rest von sexp.c als einheitliche Zufallszahl recycelt . Bisher scheint es so

• IT ist einfacher zu codieren,
• IT ist schneller und
• Sowohl IT als auch AD verlieren möglicherweise an numerischer Genauigkeit.

Ich würde mich sehr freuen, wenn jemand erklären könnte, warum R AD immer noch als einzige verfügbare Option für implementiert rexp().

Auf meinem Computer (entschuldigen Sie mein Französisch!):

> 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

Die inverse Transformation ist schlimmer. Sie sollten jedoch auf Variabilität achten. Die Einführung eines Ratenparameters führt zu einer noch größeren Variabilität der inversen Transformation:

> 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 

Hier sind die Vergleiche mit 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

Der Kilometerstand variiert also immer noch, je nach Maßstab!

Dies zitiert nur den Artikel im Abschnitt "Algorithmus LG: (Logarithmus-Methode)":

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

Es sieht also so aus, als hätten sich die Autoren für andere Methoden entschieden, um diese "Hersteller" -Einschränkung langsamer Logarithmen zu vermeiden. Vielleicht wird diese Frage dann am besten in den Stapelüberlauf verschoben, wo jemand mit Kenntnissen über die Eingeweide von R Kommentare abgeben kann.

Führen Sie dies einfach mit microbenchmark; Auf meiner Maschine ist Rs nativer Ansatz gleichmäßig schneller:

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