Transformada de Box-Muller

Transformación estadística

Visualización de la transformada de Box-Muller: los puntos coloreados en el cuadrado unitario (u ₁ , u ₂ ), dibujados como círculos, se asignan a una gaussiana 2D (z ₀ , z ₁ ), dibujada como cruces. Los gráficos en los márgenes son las funciones de distribución de probabilidad de z0 y z1. z0 y z1 no tienen límites; parecen estar en [−2,5, 2,5] debido a la elección de los puntos ilustrados. En el archivo SVG, pase el cursor sobre un punto para resaltarlo y su punto correspondiente.

La transformada de Box-Muller , de George Edward Pelham Box y Mervin Edgar Muller , ^[1] es un método de muestreo de números aleatorios para generar pares de números aleatorios independientes , estándar y distribuidos normalmente ( expectativa cero, varianza unitaria ), dada una fuente de números aleatorios distribuidos uniformemente . El método fue mencionado por primera vez explícitamente por Raymond EAC Paley y Norbert Wiener en su tratado de 1934 sobre transformadas de Fourier en el dominio complejo. ^[2] Dado el estatus de estos últimos autores y la amplia disponibilidad y uso de su tratado, es casi seguro que Box y Muller conocían bien su contenido.

La transformada de Box-Muller se expresa comúnmente en dos formas. La forma básica, tal como la dan Box y Muller, toma dos muestras de la distribución uniforme en el intervalo (0,1) y las asigna a dos muestras estándar distribuidas normalmente. La forma polar toma dos muestras de un intervalo diferente, [−1,+1] , y las asigna a dos muestras distribuidas normalmente sin el uso de funciones seno o coseno.

La transformada de Box-Muller se desarrolló como una alternativa computacionalmente más eficiente al método de muestreo de transformada inversa . ^[3] El algoritmo zigurat proporciona un método más eficiente para procesadores escalares (por ejemplo, CPU antiguas), mientras que la transformada de Box-Muller es superior para procesadores con unidades vectoriales (por ejemplo, GPU o CPU modernas). ^[4]

Forma básica

Supongamos que U ₁ y U ₂ son muestras independientes elegidas de la distribución uniforme en el intervalo unitario (0, 1) . Sea y $${\displaystyle Z_{0}=R\cos(\Theta )={\sqrt {-2\ln U_{1}}}\cos(2\pi U_{2})\,}$$ $${\displaystyle Z_{1}=R\sin(\Theta )={\sqrt {-2\ln U_{1}}}\sin(2\pi U_{2}).\,}$$

Entonces Z ₀ y Z ₁ son variables aleatorias independientes con una distribución normal estándar .

La derivación ^{[5] se basa en una propiedad de un} sistema cartesiano bidimensional , donde las coordenadas X e Y se describen mediante dos variables aleatorias independientes y distribuidas normalmente, las variables aleatorias para R ² y Θ (mostradas arriba) en las coordenadas polares correspondientes también son independientes y se pueden expresar como y $${\displaystyle R^{2}=-2\cdot \ln U_{1}\,}$$ $${\displaystyle \Theta =2\pi U_{2}.\,}$$

Como R ² es el cuadrado de la norma de la variable normal bivariada estándar ( X , Y ) , tiene una distribución de chi-cuadrado con dos grados de libertad. En el caso especial de dos grados de libertad, la distribución de chi-cuadrado coincide con la distribución exponencial y la ecuación para R ² anterior es una forma sencilla de generar la variable exponencial requerida.

Forma polar

Se utilizan dos valores distribuidos uniformemente, u y v, para producir el valor s = R ² , que también está distribuido uniformemente. Las definiciones de seno y coseno se aplican luego a la forma básica de la transformada de Box-Muller para evitar el uso de funciones trigonométricas.

La forma polar fue propuesta por primera vez por J. Bell ^[6] y luego modificada por R. Knop ^[7] . Si bien se han descrito varias versiones diferentes del método polar, aquí se describirá la versión de R. Knop porque es la más utilizada, en parte debido a su inclusión en Recetas numéricas . D. Knuth describe una forma ligeramente diferente como "Algoritmo P" en El arte de la programación informática ^[8] .

Dados u y v , independientes y uniformemente distribuidos en el intervalo cerrado [−1, +1] , establezca s = R ² = u ² + v ² . Si s = 0 o s ≥ 1 , descarte u y v , y pruebe con otro par ( u , v ) . Debido a que u y v están uniformemente distribuidos y debido a que solo se han admitido puntos dentro del círculo unitario, los valores de s también estarán uniformemente distribuidos en el intervalo abierto (0, 1) . Esto último se puede ver calculando la función de distribución acumulativa para s en el intervalo (0, 1) . Esta es el área de un círculo con radio , dividida por . A partir de esto, encontramos que la función de densidad de probabilidad tiene el valor constante 1 en el intervalo (0, 1) . De igual modo, el ángulo θ dividido por está uniformemente distribuido en el intervalo [0, 1) e independiente de s . ${\textstyle {\sqrt {s}}}$ ${\displaystyle \pi }$ ${\displaystyle 2\pi }$

Ahora identificamos el valor de s con el de U ₁ y con el de U ₂ en la forma básica. Como se muestra en la figura, los valores de y en la forma básica se pueden reemplazar con los cocientes y , respectivamente. La ventaja es que se puede evitar el cálculo directo de las funciones trigonométricas. Esto es útil cuando las funciones trigonométricas son más costosas de calcular que la división simple que reemplaza a cada una. ${\displaystyle \theta /(2\pi )}$ ${\displaystyle \cos \theta =\cos 2\pi U_{2}}$ ${\displaystyle \sin \theta =\sin 2\pi U_{2}}$ ${\displaystyle \cos \theta =u/R=u/{\sqrt {s}}}$ ${\textstyle \sin \theta =v/R=v/{\sqrt {s}}}$

Así como la forma básica produce dos desviaciones normales estándar, también lo hace este cálculo alternativo . $${\displaystyle z_{0}={\sqrt {-2\ln U_{1}}}\cos(2\pi U_{2})={\sqrt {-2\ln s}}\left({\frac {u}{\sqrt {s}}}\right)=u\cdot {\sqrt {\frac {-2\ln s}{s}}}}$$ $${\displaystyle z_{1}={\sqrt {-2\ln U_{1}}}\sin(2\pi U_{2})={\sqrt {-2\ln s}}\left({\frac {v}{\sqrt {s}}}\right)=v\cdot {\sqrt {\frac {-2\ln s}{s}}}.}$$

Contrastando las dos formas

El método polar se diferencia del método básico en que es un tipo de muestreo de rechazo . Descarta algunos números aleatorios generados, pero puede ser más rápido que el método básico porque es más sencillo de calcular (siempre que el generador de números aleatorios sea relativamente rápido) y es más robusto numéricamente. ^[9] Evitar el uso de funciones trigonométricas costosas mejora la velocidad sobre la forma básica. ^[6] Descarta 1 − π /4 ≈ 21,46% del total de pares de números aleatorios de entrada distribuidos uniformemente generados, es decir, descarta 4/ π − 1 ≈ 27,32% pares de números aleatorios distribuidos uniformemente por par de números aleatorios gaussianos generado, lo que requiere 4/ π ≈ 1,2732 números aleatorios de entrada por número aleatorio de salida.

La forma básica requiere dos multiplicaciones, 1/2 logaritmo, 1/2 raíz cuadrada y una función trigonométrica para cada variable normal. ^[10] En algunos procesadores, el coseno y el seno del mismo argumento se pueden calcular en paralelo utilizando una sola instrucción. En particular, para las máquinas basadas en Intel, se puede utilizar la instrucción de ensamblador fsincos o la instrucción expi (generalmente disponible en C como una función intrínseca ), para calcular complejos y simplemente separar las partes reales e imaginarias. $${\displaystyle \exp(iz)=e^{iz}=\cos z+i\sin z,\,}$$

Nota: Para calcular explícitamente la forma compleja-polar utilice las siguientes sustituciones en la forma general,

Dejar y luego ${\textstyle r={\sqrt {-\ln(u_{1})}}}$ ${\textstyle z=2\pi u_{2}.}$ $${\displaystyle re^{iz}={\sqrt {-\ln(u_{1})}}e^{i2\pi u_{2}}={\sqrt {-2\ln(u_{1})}}\left[\cos(2\pi u_{2})+i\sin(2\pi u_{2})\right].}$$

La forma polar requiere 3/2 multiplicaciones, 1/2 logaritmo, 1/2 raíz cuadrada y 1/2 división para cada variable normal. El efecto es reemplazar una multiplicación y una función trigonométrica con una sola división y un bucle condicional.

Truncamiento de colas

Cuando se utiliza una computadora para producir una variable aleatoria uniforme, inevitablemente tendrá algunas imprecisiones porque hay un límite inferior en lo cerca que pueden estar los números de 0. Si el generador usa 32 bits por valor de salida, el número distinto de cero más pequeño que se puede generar es . Cuando y son iguales a esto, la transformada de Box-Muller produce una desviación aleatoria normal igual a . Esto significa que el algoritmo no producirá variables aleatorias más de 6,660 desviaciones estándar de la media. Esto corresponde a una proporción de perdida debido al truncamiento, donde es la distribución normal acumulativa estándar. Con 64 bits, el límite se lleva a desviaciones estándar, para las cuales . ${\displaystyle 2^{-32}}$ ${\displaystyle U_{1}}$ ${\displaystyle U_{2}}$ ${\textstyle \delta ={\sqrt {-2\ln(2^{-32})}}\cos(2\pi 2^{-32})\approx 6.660}$ ${\displaystyle 2(1-\Phi (\delta ))\simeq 2.738\times 10^{-11}}$ ${\displaystyle \Phi (\delta )}$ ${\displaystyle \delta =9.419}$ ${\displaystyle 2(1-\Phi (\delta ))<5\times 10^{-21}}$

Implementación

C++

La transformación Box-Muller estándar genera valores de la distribución normal estándar ( es decir , desviaciones normales estándar ) con media 0 y desviación estándar 1. La implementación a continuación en C++ estándar genera valores de cualquier distribución normal con media y varianza . Si es una desviación normal estándar, entonces tendrá una distribución normal con media y desviación estándar . El generador de números aleatorios se ha inicializado para garantizar que se devuelvan nuevos valores pseudoaleatorios a partir de llamadas secuenciales a la función. ${\displaystyle \mu }$ ${\displaystyle \sigma ^{2}}$ ${\displaystyle Z}$ ${\displaystyle X=Z\sigma +\mu }$ ${\displaystyle \mu }$ ${\displaystyle \sigma }$generateGaussianNoise

#include <cmath> #include <limits> #include <random> #include <utility>    //"mu" es la media de la distribución y "sigma" es la desviación estándar. std :: pair < ​​double , double > generateGaussianNoise ( double mu , double sigma ) { constexpr double two_pi = 2.0 * M_PI ;             //inicializar el generador de números uniformes aleatorios (runif) en un rango de 0 a 1 static std :: mt19937 rng ( std :: random_device {}()); // Motor mersenne_twister_engine estándar sembrado con rd() static std :: uniform_real_distribution <> runif ( 0.0 , 1.0 );         //crea dos números aleatorios, asegúrate de que u1 sea mayor que cero double u1 , u2 ; do { u1 = runif ( rng ); } while ( u1 == 0 ); u2 = runif ( rng );                 //calcula z0 y z1 auto mag = sigma * sqrt ( -2.0 * log ( u1 )); auto z0 = mag * cos ( two_pi * u2 ) + mu ; auto z1 = mag * sin ( dos_pi * u2 ) + mu ;                             devolver std :: hacer_par ( z0 , z1 ); }  

JavaScript

función rand_normal () {   /* Sintaxis:  *  * [ x, y ] = rand_normal();  * x = rand_normal()[0];  * y = rand_normal()[1];  *  */ // Transformada de Box-Muller: sea ​​phi = 2 * Math . PI * Math . aleatorio (); sea R = Math . sqrt ( - 2 * Math . log ( Math . aleatorio () ) ); sea x = R * Math . cos ( phi ); sea y = R * Math . sin ( phi );                              // Valores de retorno: devuelve [ x , y ];    }

Julia

"""  muestra de caja (N)Genere `2N` muestras de la distribución normal estándar utilizando el método Box-Muller. """ function boxmullersample ( N ) z = Array { Float64 }( undef , N , 2 ); for i in axes ( z , 1 ) z [ i , : ] .= sincospi ( 2 * rand ()); z [ i , : ] .*= sqrt ( - 2 * log ( rand ())); end vec ( z ) end                    """  boxmullersample(n,μ,σ)Genere `n` muestras de la distribución normal con media `μ` y desviación estándar `σ` utilizando el método de Box-Muller. """ function boxmullersample ( n , μ , σ ) μ .+ σ * boxmullersample ( cld ( n , 2 ))[ 1 : n ]; fin    

Véase también

• Muestreo por transformada inversa
• Método polar de Marsaglia , transformación similar a Box-Muller, que utiliza coordenadas cartesianas, en lugar de coordenadas polares.

Referencias

1. Box, GEP; Muller, Mervin E. (1958). "Una nota sobre la generación de desviaciones normales aleatorias". Anales de estadística matemática . 29 (2): 610–611. doi : 10.1214/aoms/1177706645 . JSTOR  2237361.
2. Raymond EAC Paley y Norbert Wiener Transformadas de Fourier en el dominio complejo, Nueva York: American Mathematical Society (1934) §37.
3. Kloeden y Platen, Soluciones numéricas de ecuaciones diferenciales estocásticas , págs. 11-12
4. Howes, Lee; Thomas, David (2008). GPU Gems 3: generación y aplicación eficiente de números aleatorios con CUDA . Pearson Education, Inc. ISBN 978-0-321-51526-1.
5. Sheldon Ross, Un primer curso de probabilidad , (2002), págs. 279-281
6. Bell, James R. (1968). "Algoritmo 334: Desviaciones aleatorias normales". Comunicaciones de la ACM . 11 (7): 498. doi : 10.1145/363397.363547 .
7. Knop, R. (1969). "Observación sobre el algoritmo 334 [G5]: Desviaciones aleatorias normales". Comunicaciones de la ACM . 12 (5): 281. doi : 10.1145/362946.362996 .
8. Knuth, Donald (1998). El arte de la programación informática: Volumen 2: Algoritmos seminuméricos . pág. 122. ISBN 0-201-89684-2.
9. Everett F. Carter, Jr., La generación y aplicación de números aleatorios, Forth Dimensions (1994), vol. 16, núm. 1 y 2.
10. La evaluación de 2 π U ₁ se cuenta como una multiplicación porque el valor de 2 π se puede calcular de antemano y utilizar repetidamente.

Enlaces externos

• Weisstein, Eric W. "Transformación de Box-Muller". MathWorld .
• Cómo convertir una distribución uniforme en una distribución gaussiana (código C)