Quais são as vantagens de um gerador aleatório exponencial usando o método de Ahrens e Dieter (1972) e não por transformação inversa?

11

Minha pergunta é inspirada no gerador de números aleatórios exponenciais embutidos de R , a função rexp(). Ao tentar gerar números aleatórios distribuídos exponencialmente, muitos livros recomendam o método de transformação inversa, conforme descrito nesta página da Wikipedia . Estou ciente de que existem outros métodos para realizar essa tarefa. Em particular, o código fonte de R usa o algoritmo descrito em um artigo de Ahrens & Dieter (1972) .

Eu me convenci de que o método Ahrens-Dieter (AD) está correto. Ainda assim, não vejo o benefício de usar seu método comparado ao método de transformação inversa (TI). O AD não é apenas mais complexo de implementar que o TI. Também não parece haver um benefício na velocidade. Aqui está o meu código R para comparar os dois métodos seguidos pelos resultados.

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

Resultados:

[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

Comparando o código para os dois métodos, o AD desenha pelo menos dois números aleatórios uniformes (com a função Cunif_rand() ) para obter um número aleatório exponencial. A TI precisa apenas de um número aleatório uniforme. Presumivelmente, a equipe principal do R decidiu não implementar a TI porque supunha que assumir o logaritmo pode ser mais lento do que gerar números aleatórios mais uniformes. Entendo que a velocidade de obtenção dos logaritmos pode depender da máquina, mas pelo menos para mim o contrário é verdadeiro. Talvez haja problemas em torno da precisão numérica da TI relacionada à singularidade do logaritmo em 0? Mas então, o código fonte R sexp.crevela que a implementação do AD também perde alguma precisão numérica porque a parte a seguir do código C remove os bits iniciais do número aleatório 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 é depois reciclado como um número aleatório uniforme no restante sexp.c . Até agora, parece que

  • É mais fácil codificar a TI,
  • A TI é mais rápida e
  • o TI e o AD possivelmente perdem a precisão numérica.

Eu realmente apreciaria se alguém pudesse explicar por que o R ainda implementa o AD como a única opção disponível para rexp().

Pratyush Mais
fonte
4
Com os geradores de números aleatórios, "mais fácil de codificar" não é realmente uma consideração, a menos que você esteja fazendo isso! Velocidade e precisão são as únicas duas considerações. (Para geradores uniformes, há também o período do gerador.) Antigamente, o AD era mais rápido. Na minha caixa Linux, o AD é executado aproximadamente em 1/2 do tempo que a função invTrans, e no meu laptop em cerca de 2/3 do tempo. Você também pode usar a marca de microbench para tempos mais abrangentes.
jbowman
5
Eu sugeriria que não migrássemos. Isso me parece um tópico.
Ameba diz Reinstate Monica
1
Dado que sou incapaz de pensar em um único cenário em que rexp(n)seria o gargalo, a diferença de velocidade não é um argumento forte para a mudança (pelo menos para mim). Talvez eu esteja mais preocupado com a precisão numérica, embora não esteja claro para mim qual deles seria mais confiável numericamente.
Cliff AB
1
@amoeba Penso que "Quais seriam as vantagens de ..." seria uma reformulação que seria claramente sobre o assunto aqui e não afetaria nenhuma resposta existente. Suponho que "por que as pessoas que fizeram R decidirem fazer ..." é realmente (a) uma questão específica de software, (b) requer evidência na documentação ou telepatia, portanto, pode-se argumentar que isso é fora de tópico aqui. Pessoalmente, prefiro que a questão seja reformulada para torná-la mais clara no escopo do site, mas não vejo isso como uma razão suficientemente forte para fechá-la.
Silverfish
1
@amoeba Eu tive uma chance. Não convencido de que meu novo título sugerido seja especialmente gramatical, e talvez algumas outras partes do texto da pergunta possam mudar. Mas espero que isso seja mais claro no tópico, pelo menos, e não acho que isso invalide ou exija mudanças nas respostas.
Silverfish

Respostas:

9

No meu computador (perdoe meu francês!):

> 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

a transformação inversa faz pior. Mas você deve observar a variabilidade. A introdução de um parâmetro de taxa leva a ainda mais variabilidade para a transformação inversa:

> 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 

Aqui estão as comparações usando 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

Portanto, a milhagem ainda varia, dependendo da balança!

Xi'an
fonte
2
No meu laptop, os horários são tão parecidos com os do OP que suspeito que tenhamos a mesma máquina (ou pelo menos o mesmo processador). Mas acho que seu ponto aqui é que a vantagem de velocidade observada depende da plataforma e, dada a diferença mínima, não há uma vantagem clara em relação à velocidade.
Cliff AB
4
Você poderia, talvez, executar um microbenchmark?
Firebug
2
rexp-log(runif())5.27±0.02Rlogrunif
7

Isso está apenas citando o artigo na seção "Algorithm LG: (Logarithm method)":

X=ALOG(REGOL(IR))μμμu

Portanto, parece que os autores optaram por outros métodos para evitar essa limitação do "fabricante" dos logaritmos lentos. Talvez essa questão seja melhor movida para o stackoverflow, onde alguém com conhecimento sobre as intenções de R pode comentar.

Alex R.
fonte
6

Apenas executando isso com microbenchmark; na minha máquina, a abordagem nativa de R é uniformemente mais rápida:

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

insira a descrição da imagem aqui

MichaelChirico
fonte