¿Ventajas de Box-Muller sobre método inverso de la FCD para simular la distribución Normal?

Question

Con el fin de simular una distribución normal a partir de un conjunto uniforme de las variables, hay varias técnicas:

1. El Box-Muller algoritmo, en el que dos muestras independientes uniforme variables en $(0,1)$, y los transforma en dos independientes estándar de las distribuciones normales a través de: $$ Z_0 = \sqrt{-2\text{ln}U_1}\text{cos}(2\pi U_0)\\ Z_1 = \sqrt{-2\text{ln}U_1}\text{pecado}(2\pi U_0) $$

2. la CDF método, donde uno se puede equiparar el normal cdf $(F(Z))$ a un Uniforme de la variable aleatoria: $$ F(Z) = U $$ y derivar $$Z = F^{-1}(U)$$

Mi pregunta es: cual es computacionalmente más eficiente? Yo creo que el último método, pero la mayoría de los artículos que he leído uso de Box-Muller - por qué?

Información adicional:

La inversa de la normal CDF es saber y dada por:

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

Por lo tanto: $$ Z = F^{-1}(U)\; =\; \sqrt2\;\operatorname{fer}^{-1}(2U - 1), \quad U\en(0,1).$$

Answer 1

Desde una perspectiva puramente probabilístico perspectiva, ambos enfoques son correctos y por lo tanto equivalente. Desde una perspectiva algorítmica, la comparación debe considerar tanto la precisión y la computación costo.

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

• generar $$U_1,U_2\stackrel{\text{iid}}{\sim}\mathcal{U}(-1,1)$$ until $S=U_1^2+U_2^2\le 1$
• tome $Z=\sqrt{-2\log(S)/S}$ y definen $$X_1=ZU_1\,,\ X_2=Z U_2$$

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

A seguir en los comentarios realizados por whuber, la comparación de rnorm(N)y qnorm(runif(N))está en la ventaja de la inversa de la cdf, 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:

Después de un comentario de Radford a Neal en mi blog, quiero señalar que el valor por defecto rnorm en R hace uso de la inversión de método, de ahí que la comparación anterior se refleja en la interfaz y no en el método de simulación de sí mismo! Para citar la R de documentación sobre el 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.