Zalety Box-Mullera nad odwrotną metodą CDF do symulacji rozkładu normalnego?

Aby zasymulować rozkład normalny z zestawu zmiennych jednorodnych, istnieje kilka technik:

Algorytm Boxa-Mullera , w którym jeden próbkuje dwa niezależne jednolite zmienia się na $(0,1)$ i przekształca je w dwa niezależne standardowe rozkłady normalne poprzez:
$Z_{0} = \sqrt{- 2 ln U_{1}} cos (2 π U_{0}) Z_{1} = \sqrt{- 2 ln U_{1}} sin (2 π U_{0})$ $Z_0 = \sqrt{-2\text{ln}U_1}\text{cos}(2\pi U_0)\\ Z_1 = \sqrt{-2\text{ln}U_1}\text{sin}(2\pi U_0)$
metoda CDF , w której można zrównać normalny cdf $(F(Z))$ do wariantu Uniform:
$F (Z) = U$ $F(Z) = U$ i wyprowadzić $Z = F^{- 1} (U)$ $Z = F^{-1}(U)$

Moje pytanie brzmi: która jest bardziej wydajna obliczeniowo? Myślałem, że to ta druga metoda - ale większość artykułów, które czytam, używa Box-Mullera - dlaczego?

Dodatkowe informacje:

Odwrotność normalnego CDF jest znana i podana przez:

F^{- 1} (Z) = \sqrt{2} \erf^{- 1} (2 Z - 1), Z \in (0, 1) .

$F^{-1}(Z)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2Z - 1), \quad Z\in(0,1).$

Stąd:

Z = F^{- 1} (U) = \sqrt{2} \erf^{- 1} (2 U - 1), U \in (0, 1) .

$Z = F^{-1}(U)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2U - 1), \quad U\in(0,1).$

normal-distribution simulation uniform

— użytkownik2350366
źródło

Jaka jest odwrotność normalnego cdf? Nie można tego obliczyć analitycznie, tylko jeśli oryginalny CDF jest aproksymowany za pomocą częściowej funkcji liniowej.

— Artem Sobolev

Czy te dwa nie są ze sobą blisko powiązane? Uważam, że Box Muller to szczególny przypadek generacji 2-zmiennej.

— ttnphns

Cześć Barmaley, powyżej dodałem więcej informacji. Inverse CDF ma wyrażenie - jednak

musi być obliczony obliczeniowo - więc może dlatego preferowane jest pole Muller? Zakładałem, że

będzie obliczany w tabelach odnośników, podobnie jak wartości

? Czy zatem nie jest to o wiele bardziej kosztowne obliczeniowo? Mogę się jednak mylić.

{erf}^{- 1}

$\text{erf}^{-1}$

{erf}^{1}

$\text{erf}^1$

sin

$\text{sin}$

cosine

$\text{cosine}$

— user2350366,

Istnieją wersje Box-Mullera bez grzechu i cosinusa.

— Xi'an

@Dilip W przypadku aplikacji o bardzo niskiej precyzji, takich jak grafika komputerowa, sinus i cosinus można rzeczywiście zoptymalizować za pomocą odpowiednich tabel odnośników. Jednak w przypadku zastosowań statystycznych taka optymalizacja nigdy nie jest wykorzystywana. Ostatecznie nie jest tak naprawdę trudniej obliczyć

niż

lub

, ale w nowoczesnych systemach komputerowych elementarne funkcje związane z

tym funkcje triggera - mają tendencję do optymalizacji (

były podstawowymi instrukcjami dawno temu na Intelie 8087 chip!), Podczas gdy erf jest niedostępny lub został zakodowany na wyższym (= wolniejszym) poziomie.

{erf}^{- 1}

$\text{erf}^{-1}$

\log

$\log$

sqrt

$\text{sqrt}$

\exp

$\exp$

\cos

$\cos$

\log

$\log$

— whuber

Z czysto probabilistycznego punktu widzenia oba podejścia są poprawne, a zatem równoważne. Z punktu widzenia algorytmu porównanie musi uwzględniać zarówno precyzję, jak i koszt obliczeniowy.

Box-Muller opiera się na jednolitym generatorze i kosztuje mniej więcej tyle samo, co ten jednolity generator. Jak wspomniano w moim komentarzu, możesz uciec bez połączeń sinusoidalnych lub cosinusowych, jeśli nie bez logarytmu:

generować aż $U_{1}, U_{2} \overset{iid}{\sim} U (- 1, 1)$ $U_1,U_2\stackrel{\text{iid}}{\sim}\mathcal{U}(-1,1)$ $S=U_1^2+U_2^2\le 1$
weź i zdefiniuj $Z=\sqrt{-2\log(S)/S}$ $X_{1} = Z U_{1}, X_{2} = Z U_{2}$ $X_1=ZU_1\,,\ X_2=Z U_2$

Ogólny algorytm inwersji wymaga na przykład wywołania odwrotnego normalnego pliku cdf qnorm(runif(N)) w R, co może być bardziej kosztowne niż powyższe i, co ważniejsze, może zawieść w ogonach pod względem precyzji, chyba że funkcja kwantylu jest dobrze zakodowana.

Aby podążać za komentarzami Whubera , porównanie rnorm(N)i qnorm(runif(N))przewaga odwrotnego cdf, zarówno w czasie wykonania:

> system.time(qnorm(runif(10^8)))
sutilisateur     système      écoulé
 10.137           0.120      10.251 
> system.time(rnorm(10^8))
utilisateur     système      écoulé
 13.417           0.060      13.472` `

i pod względem dopasowania do ogona: enter image description here

Po komentarzu Radforda Neala na moim blogu , chcę zauważyć, że domyślna wartość rnormw R wykorzystuje metodę inwersji, stąd powyższe porównanie odzwierciedla interfejs, a nie samą metodę symulacji! Aby zacytować dokumentację R dotyczącą RNG:

‘normal.kind’ can be ‘"Kinderman-Ramage"’, ‘"Buggy
 Kinderman-Ramage"’ (not for ‘set.seed’), ‘"Ahrens-Dieter"’,
 ‘"Box-Muller"’, ‘"Inversion"’ (the default), or ‘"user-supplied"’.
 (For inversion, see the reference in ‘qnorm’.)  The
 Kinderman-Ramage generator used in versions prior to 1.7.1 (now
 called ‘"Buggy"’) had several approximation errors and should only
 be used for reproduction of old results.  The ‘"Box-Muller"’
 generator is stateful as pairs of normals are generated and
 returned sequentially.  The state is reset whenever it is selected
 (even if it is the current normal generator) and when ‘kind’ is
 changed.

— Xi'an
źródło

(1) Jak są

\log

$\log$

\sqrt{}

$\sqrt{}$

Φ^{- 1}

$\Phi^{-1}$

Φ^{- 1}

$\Phi^{-1}$

X_{1}

$X_1$

X_{2}

$X_2$

U_{i}

$U_i$

- 1

$-1$

1

$1$

0

$0$

1

$1$

R 3.0.2rowSums

S

$S$ qnorm(runif(N))InverseCDF[NormalDistribution[], #] &

Zgadzam się, qnorm(runif(N))jest nawet o 20% szybszy niżrnorm(N)

— Xi'an

Φ^{- 1}

$\Phi^{-1}$

\sin

$\sin$

\cos

$\cos$

Dla porównania za pomocą i7-3740QM @ 2,7 Ghz i R 3.12 dla następujących wywołań:

RNGkind(kind = NULL, normal.kind = 'Inversion');At <- microbenchmark(A <- rnorm(1e5, 0, 1), times = 100L);RNGkind(kind = NULL, normal.kind = 'Box-Muller');Bt <- microbenchmark(B <- rnorm(1e5, 0, 1), times = 100L)

Dostaję mean 11.38363 median 11.18718inwersję i mean 13.00401 median 12.48802dla Box-Mullera

— Avraham