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

11

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

Pratyush Mehr
quelle
4
Bei Zufallszahlengeneratoren spielt "einfacher zu codieren" keine Rolle, es sei denn, Sie sind derjenige, der dies tut! Geschwindigkeit und Genauigkeit sind die einzigen beiden Überlegungen. (Bei einheitlichen Generatoren gibt es auch die Periode des Generators.) Früher war AD schneller. Auf meiner Linux-Box läuft AD in ungefähr der Hälfte der Zeit, die Ihre invTrans-Funktion ausführt, und auf meinem Laptop in ungefähr 2/3 der Zeit. Möglicherweise möchten Sie Microbenchmark auch für umfassendere Timings verwenden.
Jbowman
5
Ich würde vorschlagen, dass wir es nicht migrieren. Dies scheint mir ein Thema zu sein.
Amöbe sagt Reinstate Monica
1
Da ich mir kein einziges Szenario ausdenken kann, in dem rexp(n)es zu einem Engpass kommen würde, ist der Geschwindigkeitsunterschied (zumindest für mich) kein starkes Argument für Veränderungen. Ich bin möglicherweise mehr besorgt über die numerische Genauigkeit, obwohl mir nicht klar ist, welche numerisch zuverlässiger wäre.
Cliff AB
1
@amoeba Ich denke, dass "Was wären die Vorteile von ..." eine Umformulierung wäre, die hier eindeutig zum Thema gehört und keine bestehenden Antworten beeinflusst. Ich nehme an, "Warum haben sich die Leute, die R dazu gebracht haben, zu tun ..." ist wirklich (a) eine softwarespezifische Frage, (b) erfordert entweder Beweise in der Dokumentation oder Telepathie, könnte also hier wohl nicht zum Thema gehören. Persönlich wäre es mir lieber, wenn die Frage umformuliert würde, um sie im Rahmen der Website klarer zu machen, aber ich sehe dies nicht als ausreichend starken Grund, sie zu schließen.
Silverfish
1
@ Amöbe Ich habe es versucht. Ich bin nicht überzeugt, dass mein vorgeschlagener neuer Titel besonders grammatikalisch ist, und vielleicht könnten einige andere Teile des Fragentextes eine Änderung vertragen. Aber ich hoffe, dass dies zumindest klarer zum Thema gehört, und ich denke nicht, dass es ungültig macht oder Änderungen an beiden Antworten erfordert.
Silverfish

Antworten:

9

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!

Xi'an
quelle
2
Auf meinem Laptop stimmen die Zeiten so genau mit den OPs überein, dass ich vermute, dass wir dieselbe Maschine (oder zumindest denselben Prozessor) haben. Aber ich denke, Ihr Punkt hier ist, dass der beobachtete Geschwindigkeitsvorteil plattformabhängig ist, und angesichts des minimalen Unterschieds gibt es keinen klaren Vorteil gegenüber dem anderen in Bezug auf die Geschwindigkeit.
Cliff AB
4
Könnten Sie vielleicht microbenchmarkstattdessen eine durchführen?
Firebug
2
rexp-log(runif())5.27±0,02Rlogrunif
7

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

X.=- -EINL.ÖG(R.E.GÖL.(ichR.))μμμ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.

Alex R.
quelle
6

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

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

Geben Sie hier die Bildbeschreibung ein

MichaelChirico
quelle