Mi pregunta está inspirada en el generador de números aleatorios exponencial incorporado de R , la función rexp(). Al intentar generar números aleatorios distribuidos exponencialmente, muchos libros de texto recomiendan el método de transformación inversa como se describe en esta página de Wikipedia . Soy consciente de que hay otros métodos para realizar esta tarea. En particular, el código fuente de R utiliza el algoritmo descrito en un artículo de Ahrens y Dieter (1972) .
Me he convencido de que el método Ahrens-Dieter (AD) es correcto. Aún así, no veo el beneficio de usar su método en comparación con el método de transformación inversa (IT). AD no solo es más complejo de implementar que TI. Tampoco parece haber un beneficio de velocidad. Aquí está mi código R para comparar ambos métodos seguidos de los 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
Al comparar el código de los dos métodos, AD dibuja al menos dos números aleatorios uniformes (con la función Cunif_rand() ) para obtener un número aleatorio exponencial. Solo necesita un número aleatorio uniforme. Presumiblemente, el equipo central de R decidió no implementar TI porque asumió que tomar el logaritmo puede ser más lento que generar números aleatorios más uniformes. Entiendo que la velocidad de tomar logaritmos puede depender de la máquina, pero al menos para mí lo contrario es cierto. ¿Quizás haya problemas en torno a la precisión numérica de TI que tiene que ver con la singularidad del logaritmo en 0? Pero entonces, el
código fuente R sexp.crevela que la implementación de AD también pierde algo de precisión numérica porque la siguiente parte del código C elimina los bits iniciales del número aleatorio 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 más tarde se recicla como un número aleatorio uniforme en el resto de sexp.c . Hasta ahora, parece que
- Es más fácil codificar,
- Es más rápido y
- tanto IT como AD posiblemente pierdan precisión numérica.
Realmente agradecería si alguien pudiera explicar por qué R todavía implementa AD como la única opción disponible para rexp().
fuente
rexp(n)sería el cuello de botella, la diferencia de velocidad no es un argumento fuerte para el cambio (al menos para mí). Podría estar más preocupado por la precisión numérica, aunque no me queda claro cuál sería más confiable numéricamente.Respuestas:
En mi computadora (¡perdón por mi francés!):
la transformación inversa empeora. Pero debes tener cuidado con la variabilidad. La introducción de un parámetro de velocidad conduce a una variabilidad aún mayor para la transformación inversa:
Aquí están las comparaciones usando
rbenchmark:¡Entonces el kilometraje todavía varía, dependiendo de la escala!
fuente
microbenchmarken su lugar?rexp, 3% tiempo de usuario más corto para-log(runif()), y sin variabilidad con un parámetro de velocidad ( segundos totales). Todos estamos asumiendo implícitamente que está logrando tiempos y es comparable a lo que uno obtendría con una subrutina Fortran.RlogrunifEsto es solo citando el artículo en la sección "Algoritmo LG: (Método de logaritmo)":
Por lo tanto, parece que los autores optaron por otros métodos para evitar esta limitación "fabricante" de logaritmos lentos. Quizás esta pregunta se traslade mejor a stackoverflow donde alguien con conocimiento en las entrañas de R pueda comentar.
fuente
Solo ejecuta esto con
microbenchmark; en mi máquina, el enfoque nativo de R es uniformemente más rápido:Por el bien de la novedad, aquí nos aseguramos de que no se deba totalmente a tener :λ=1
fuente