Quels sont les avantages d'un générateur exponentiel aléatoire utilisant la méthode d'Ahrens et Dieter (1972) plutôt que par transformée inverse?

Ma question est inspirée du générateur exponentiel de nombres aléatoires intégré de R , la fonction rexp(). Lorsque vous essayez de générer des nombres aléatoires distribués de façon exponentielle, de nombreux manuels recommandent la méthode de transformation inverse, comme indiqué dans cette page Wikipedia . Je suis conscient qu'il existe d'autres méthodes pour accomplir cette tâche. En particulier, le code source de R utilise l'algorithme décrit dans un article d'Ahrens & Dieter (1972) .

Je me suis convaincu que la méthode Ahrens-Dieter (AD) est correcte. Pourtant, je ne vois pas l'avantage d'utiliser leur méthode par rapport à la méthode de transformation inverse (IT). AD n'est pas seulement plus complexe à mettre en œuvre que l'informatique. Il ne semble pas non plus y avoir de gain de vitesse. Voici mon code R pour comparer les deux méthodes suivies des résultats.

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

Résultats:

[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

En comparant le code pour les deux méthodes, AD dessine au moins deux nombres aléatoires uniformes (avec la fonction Cunif_rand() ) pour obtenir un nombre aléatoire exponentiel. Le service informatique n'a besoin que d'un nombre aléatoire uniforme. On peut supposer que l' équipe principale R a décidé de ne pas mettre en œuvre l'informatique, car elle a supposé que la prise du logarithme peut être plus lente que la génération de nombres aléatoires plus uniformes. Je comprends que la vitesse de prise des logarithmes peut être dépendante de la machine, mais au moins pour moi, le contraire est vrai. Peut-être y a-t-il des problèmes de précision numérique de l'informatique liés à la singularité du logarithme à 0? Mais alors, le code source R sexp.crévèle que la mise en œuvre de AD perd également une certaine précision numérique car la partie suivante du code C supprime les bits de tête du nombre aléatoire uniforme 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 est ensuite recyclé comme un nombre aléatoire uniforme dans la suite de sexp.c . Jusqu'à présent, il semble que

• IL est plus facile de coder,
• L'informatique est plus rapide et
• IT et AD perdent probablement la précision numérique.

J'apprécierais vraiment si quelqu'un pouvait expliquer pourquoi R implémente toujours AD comme seule option disponible pour rexp().

Sur mon ordinateur (pardonnez mon français!):

> 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

la transformée inverse fait pire. Mais vous devez faire attention à la variabilité. L'introduction d'un paramètre de taux entraîne encore plus de variabilité pour la transformée inverse:

> 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 

Voici les comparaisons utilisant 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

Le kilométrage varie donc toujours selon l'échelle!

Ceci ne fait que citer l'article dans la section "Algorithme LG: (méthode du logarithme)":

Dans FORTRAN, l'algorithme est mieux programmé directement comme évitant tout appel de sous-programme. La performance était de 504 sec par échantillon. De ce temps, 361 sec ont été absorbés par la routine de logarithme du fabricant et 105 sec par le générateur pour REGOL de la variable uniformément répartie . Il était donc inutile d'essayer d'améliorer la vitesse en écrivant une fonction assembleur qui devrait utiliser les deux mêmes sous-programmes.μ μ μ $X=-ALOG(REGOL(IR))$$\mu$$\mu$$\mu$$u$

Il semble donc que les auteurs aient opté pour d'autres méthodes pour éviter cette limitation «constructeur» des logarithmes lents. Peut-être que cette question est alors mieux déplacée vers stackoverflow où quelqu'un ayant des connaissances sur les tripes de R peut commenter.

Il suffit de lancer cela avec microbenchmark; sur ma machine, l'approche native de R est uniformément plus rapide:

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

Par souci de nouveauté, assurez-vous que ce n'est pas entièrement dû à la présence de :$\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()