Moje pytanie jest inspirowane wbudowanym generatorem wykładniczej liczby losowej R. , funkcją rexp(). Podczas próby generowania wykładniczych liczb losowych rozkładanych wykładniczo wiele podręczników zaleca metodę transformacji odwrotnej opisaną na tej stronie Wikipedii . Wiem, że istnieją inne metody realizacji tego zadania. W szczególności kod źródłowy R korzysta z algorytmu przedstawionego w artykule Ahrensa i Dietera (1972) .
Przekonałem się, że metoda Ahrensa-Dietera (AD) jest poprawna. Mimo to nie widzę korzyści z zastosowania ich metody w porównaniu z metodą transformacji odwrotnej (IT). AD jest nie tylko bardziej skomplikowany do wdrożenia niż IT. Wydaje się, że nie ma również korzyści z prędkości. Oto mój kod R do porównania obu metod, a następnie wyników.
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)))
Wyniki:
[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
Porównując kod dla dwóch metod, AD rysuje co najmniej dwie jednolite liczby losowe (z funkcją Cunif_rand() ), aby uzyskać jedną wykładniczą liczbę losową. IT potrzebuje tylko jednej jednolitej liczby losowej. Prawdopodobnie główny zespół R zdecydował się nie wdrażać technologii informatycznych, ponieważ zakładał, że przyjęcie logarytmu może być wolniejsze niż generowanie bardziej jednolitych liczb losowych. Rozumiem, że szybkość przyjmowania logarytmów może zależeć od maszyny, ale przynajmniej dla mnie jest odwrotnie. Być może istnieją problemy związane z dokładnością liczbową IT związaną z osobliwością logarytmu przy 0? Ale potem
kod źródłowy R sexp.cujawnia, że implementacja AD również traci pewną dokładność liczbową, ponieważ następna część kodu C usuwa wiodące bity z jednolitej liczby losowej 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 jest następnie zawracana jako jednolity losowej liczby w pozostałej sexp.c . Jak dotąd wydaje się, jakby
- Łatwiej jest kodować,
- IT jest szybsze i
- zarówno IT, jak i AD prawdopodobnie tracą dokładność liczbową.
Byłbym naprawdę wdzięczny, gdyby ktoś mógł wyjaśnić, dlaczego R nadal implementuje AD jako jedyną dostępną opcję dla rexp().
rexp(n)byłoby wąskie gardło, różnica prędkości nie jest silnym argumentem za zmianą (przynajmniej dla mnie). Być może bardziej martwię się o dokładność numeryczną, chociaż nie jest dla mnie jasne, który z nich byłby bardziej wiarygodny numerycznie.
