Ventajas de Box-Muller sobre el método CDF inverso para simular la distribución normal?

Para simular una distribución normal a partir de un conjunto de variables uniformes, existen varias técnicas:

El algoritmo Box-Muller , en el que se toman muestras de dos variables uniformes independientes en $(0,1)$ y se transforman en dos distribuciones normales estándar independientes a través de:
$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)$
El método CDF , donde se puede equiparar el cdf normal a una variante uniforme: y derivar $(F(Z))$
$F (Z) = U$ $F(Z) = U$ $Z = F^{- 1} (U)$ $Z = F^{-1}(U)$

Mi pregunta es: ¿cuál es computacionalmente más eficiente? Creo que es el último método, pero la mayoría de los documentos que leo usan Box-Muller, ¿por qué?

Información Adicional:

El inverso del CDF normal es conocido y dado por:

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

Por lo tanto:

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 usuario2350366
fuente

¿Cuál es el inverso del cdf normal? No se puede calcular analíticamente, solo si el CDF original se aproxima con una función lineal por partes.

Artem Sobolev

¿No están los dos estrechamente relacionados? Box Muller, creo, es un caso particular para la generación de 2 variantes.

ttnphns

Hola Barmaley, he agregado más información arriba. El CDF inverso tiene una expresión, sin embargo, el

debe calcularse computacionalmente, por lo que podría ser por eso que se prefiere el cuadro Muller. Supuse que

se calcularía en tablas de búsqueda, al igual que los valores de

. Por lo tanto, ¿no es mucho más costoso computacionalmente? Sin embargo, puedo estar equivocado.

{erf}^{- 1}

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

{erf}^{1}

$\text{erf}^1$

sin

$\text{sin}$

cosine

$\text{cosine}$

user2350366

Hay versiones de Box-Muller sin pecado y coseno.

Xi'an

@Dilip Para aplicaciones de muy baja precisión, como gráficos de computadora, el seno y el coseno pueden optimizarse mediante el uso de tablas de búsqueda adecuadas. Sin embargo, para aplicaciones estadísticas, dicha optimización nunca se usa. En última instancia, no es realmente más difícil calcular

que

, pero en los sistemas informáticos modernos, las funciones elementales relacionadas con

--incluidas las funciones trigonométricas-- tienden a optimizarse (

fueron instrucciones básicas en Intel) 8087 chip!), Mientras que erf no está disponible o se ha codificado a un nivel más alto (= más lento).

{erf}^{- 1}

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

\log

$\log$

sqrt

$\text{sqrt}$

\exp

$\exp$

\cos

$\cos$

\log

$\log$

whuber

Respuestas:

Desde una perspectiva puramente probabilística, ambos enfoques son correctos y, por lo tanto, equivalentes. Desde una perspectiva algorítmica, la comparación debe considerar tanto la precisión como el costo informático.

Box-Muller depende de un generador uniforme y cuesta casi lo mismo que este generador uniforme. Como se mencionó en mi comentario, puede escapar sin llamadas de seno o coseno, si no sin el logaritmo:

generar hasta $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$
tomar y define $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$

El algoritmo de inversión genérico requiere la llamada al cdf normal inverso, por ejemplo qnorm(runif(N))en R, que puede ser más costoso que el anterior y, lo que es más importante, puede fallar en las colas en términos de precisión, a menos que la función cuantil esté bien codificada.

Para seguir los comentarios hechos por whuber , la comparación de rnorm(N)y qnorm(runif(N))es una ventaja del cdf inverso, tanto en tiempo de ejecución:

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

y en términos de ajuste en la cola: enter image description here

Después de un comentario de Radford Neal en mi blog , quiero señalar que el valor predeterminado rnormen R utiliza el método de inversión, por lo tanto, la comparación anterior se refleja en la interfaz y no en el método de simulación en sí. Para citar la documentación de R en 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
fuente

\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$

whuber

R 3.0.2rowSums

S

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

Estoy de acuerdo, qnorm(runif(N))es incluso un 20% más rápido quernorm(N)

Xi'an

Φ^{- 1}

$\Phi^{-1}$

\sin

$\sin$

\cos

$\cos$ sería tan eficiente, también los usaría en lugar del muestreo de rechazo.

whuber

A modo de comparación, utilizando un i7-3740QM @ 2.7Ghz y R 3.12, para las siguientes llamadas:

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)

obtengo mean 11.38363 median 11.18718por inversión y mean 13.00401 median 12.48802por Box-Muller

Avraham el