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()
.
fonte
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.Respostas:
No meu computador (perdoe meu francês!):
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:
Aqui estão as comparações usando
rbenchmark
:Portanto, a milhagem ainda varia, dependendo da balança!
fonte
microbenchmark
?rexp
-log(runif())
R
log
runif
Isso está apenas citando o artigo na seção "Algorithm LG: (Logarithm method)":
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.
fonte
Apenas executando isso com
microbenchmark
; na minha máquina, a abordagem nativa de R é uniformemente mais rápida:fonte