Transformada de Box-Muller

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, normalmente distribuidos ( expectativa cero , varianza unitaria ), dada una fuente de números aleatorios distribuidos . El método fue mencionado explícitamente por primera vez por Raymond EAC Paley y Norbert Wiener en su tratado de 1934 sobre las 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 de dos formas. La forma básica dada por 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 por 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 1$ y $U 2$ son muestras independientes elegidas de la distribución uniforme en el intervalo unitario $(0, 1)$ . Dejar

Z_{0}=R\cos(\Theta )={\sqrt {-2\ln U_{1}}}\cos(2\pi U_{2})\,

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 normalmente distribuidas, las variables aleatorias para $R 2$ y $Θ$ (mostradas arriba) en la polar correspondiente. Las coordenadas también son independientes y se pueden expresar como

R^{2}=-2\cdot \ln U_{1}\,

\Theta =2\pi U_{2}.\,

Debido a que $R 2$ es el cuadrado de la norma de la variable normal bivariada estándar $(X, Y)$ , tiene la distribución chi-cuadrado con dos grados de libertad. En el caso especial de dos grados de libertad, la distribución chi-cuadrado coincide con la distribución exponencial , y la ecuación para $R 2$ 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 2$ , que también está distribuido uniformemente. Luego, las definiciones de seno y coseno se aplican 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, la versión de R. Knop se describirá aquí porque es la más utilizada, en parte debido a su inclusión en Numerical Recipes . D. Knuth describe una forma ligeramente diferente como "Algoritmo P" en The Art of Computer Programming . ^[8]

Dados $u$ y $v$ , independientes y distribuidos uniformemente en el intervalo cerrado $[-1, +1]$ , establezca $s = R 2 = u 2 + v 2$ . Si $s = 0$ o $s \geq 1$ , descarte u y v , y pruebe con otro par $(u, v)$ . Debido a que $u$ y $v$ están distribuidos uniformemente y debido a que sólo se han admitido puntos dentro del círculo unitario, los valores de s también estarán distribuidos uniformemente 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)$ . Asimismo, el ángulo θ dividido por está distribuido uniformemente en el intervalo $[0, 1)$ e independiente de $s$ . ${\estilo de texto {\sqrt {s}}}$ $\pi$ $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 las razones 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 única que reemplaza a cada una. $\theta /(2\pi )$ $\cos \theta =\cos 2\pi U_ {2}$ $\sin \theta =\sin 2\pi U_ {2}$ $\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.

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

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 numéricamente más robusto. ^[9] Evitar el uso de costosas funciones trigonométricas mejora la velocidad con respecto a la forma básica. ^[6] Descarta $1 - π /4 \approx 21,46%$ del total de pares de números aleatorios distribuidos uniformemente generados, es decir, descarta $4/ π - 1 \approx 27,32%$ de pares de números aleatorios distribuidos uniformemente por cada par de números aleatorios gaussianos generados, lo que requiere $4/ π \approx 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 máquinas basadas en Intel, se puede usar la instrucción de ensamblador fsincos o la instrucción expi (generalmente disponible en C como una función intrínseca ), para calcular complejos

\exp(iz)=e^{iz}=\cos z+i\sin z,\,

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

deja y luego ${\textstyle r={\sqrt {-\ln(u_{1})}}}$ ${\textstyle z=2\pi u_{2}.}$

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 existe un límite inferior sobre qué tan cerca pueden estar los números de 0. Si el generador usa 32 bits por valor de salida, el número más pequeño distinto de cero que puede ser generado 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 a más de 6,660 desviaciones estándar de la media. Esto corresponde a una proporción de la pérdida debido al truncamiento, donde es la distribución normal acumulativa estándar. Con 64 bits el límite se lleva a las desviaciones estándar, para las cuales . $2^{-32}$ ${\ Displaystyle U_ {1}}$ ${\ Displaystyle U_ {2}}$ ${\textstyle \delta ={\sqrt {-2\ln(2^{-32})}}\cos(2\pi 2^{-32})\aproximadamente 6,660}$ $2(1-\Phi (\delta ))\simeq 2.738\times 10^{-11}$ $\Phi (\delta)$ $\delta =9.419$ $2(1-\Phi (\delta ))<5\times 10^{-21}$

Implementación

C++

La transformada estándar de Box-Muller genera valores a partir de la distribución normal estándar ( es decir, desviaciones normales estándar ) con media 0 y desviación estándar 1 . La siguiente implementació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 diseñado para garantizar que se devuelvan nuevos valores pseudoaleatorios a partir de llamadas secuenciales a la función. $\mu$ $\sigma ^{2}$ $Z$ $X=Z\sigma +\mu$ $\mu$ $\sigma$ generateGaussianNoise

#include <cmath> #include <límites> #include <aleatorio> #include <utilidad>    //"mu" es la media de la distribución y "sigma" es la desviación estándar. std :: par < doble , doble > generar ruido gaussiano ( doble mu , doble sigma ) { constexpr doble dos_pi = 2.0 * M_PI ;             //inicializa el generador de números uniformes aleatorios (runif) en un rango de 0 a 1 static std :: mt19937 rng ( std :: random_device {}()); // 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 ; hacer { u1 = runif ( rng ); } mientras ( 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 ;                             return std :: make_pair ( z0 , z1 ); }

Julia

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

Ver también

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

Referencias

^ Caja, GEP; Müller, Mervin E. (1958). "Una nota sobre la generación de desviaciones normales aleatorias". Los anales de la estadística matemática . 29 (2): 610–611. doi : 10.1214/aoms/1177706645 . JSTOR 2237361.
^ Raymond EAC Paley y Norbert Wiener Fourier transforman en el dominio complejo, Nueva York: American Mathematical Society (1934) §37.
^ Kloeden y Platen, Soluciones numéricas de ecuaciones diferenciales estocásticas , págs. 11-12
^ Howes, Lee; Tomás, David (2008). GPU Gems 3: generación y aplicación eficiente de números aleatorios mediante CUDA . Pearson Educación, Inc. ISBN 978-0-321-51526-1.
^ Sheldon Ross, Un primer curso de probabilidad , (2002), págs. 279–281
^ ab Bell, James R. (1968). "Algoritmo 334: desviaciones aleatorias normales". Comunicaciones de la ACM . 11 (7): 498. doi : 10.1145/363397.363547 .
^ 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 .
^ Knuth, Donald (1998). El arte de la programación informática: Volumen 2: Algoritmos seminuméricos . pag. 122.ISBN 0-201-89684-2.
^ Everett F. Carter, Jr., La generación y aplicación de números aleatorios, Forth Dimensions (1994), vol. 16, núms. 1 y 2.
^ La evaluación de 2 $π$ U ₁ se cuenta como una multiplicación porque el valor de 2 $π$ se puede calcular con anticipación y usarse repetidamente.

enlaces externos

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