Transformación de Box-Muller

En matemáticas , la transformada de Box-Muller , introducida por 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 con distribución normal ( esperanza cero , varianza unitaria ), a partir de una fuente de números aleatorios con distribución uniforme . El método fue mencionado explícitamente por primera vez por Raymond EAC Paley y Norbert Wiener en su tratado de 1934 sobre transformadas de Fourier en el dominio complejo. ^{[ 2 ]} Dado el prestigio 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${\displaystyle (0,1)}$y los asigna a dos muestras estándar con distribución normal. La forma polar toma dos muestras de un intervalo diferente,${\displaystyle [-1,1]}$y las asigna a dos muestras con distribución normal sin utilizar funciones seno o coseno.

La transformada de Box-Muller se desarrolló como una alternativa computacionalmente más eficiente al método de muestreo de la transformada inversa . ^{[ 3 ]} El algoritmo de 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 ]}

Suponer${\displaystyle U_{1}}$y${\displaystyle U_{2}}$son muestras independientes elegidas de la distribución uniforme en el intervalo unitario${\displaystyle (0,1)}$. Colocar

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

y

${\displaystyle Z_{1}={\sqrt {-2\ln U_{1}}}\sin(2\pi U_{2})=R\sin(\Theta ).}$

Entonces${\displaystyle Z_{0}}$y${\displaystyle Z_{1}}$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 ${\displaystyle X}$y${\displaystyle Y}$Las coordenadas se describen mediante dos variables aleatorias independientes y con distribución normal; las variables aleatorias para R ² y Θ (mostradas arriba) en las coordenadas polares correspondientes también son independientes y pueden expresarse como $${\displaystyle R^{2}=-2\cdot \ln U_{1}\,}$$ y $${\displaystyle \Theta =2\pi U_{2}.\,}$$

Dado que R² es el ^cuadrado de la norma de la variable normal bivariada estándar ( X , Y ) , tiene una 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² ^anterior es una forma sencilla de generar la variable exponencial requerida.

La forma polar fue propuesta inicialmente por James Bell ^{[ 6 ]} y posteriormente 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, ya que 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 ] , definimos s = R ² = u ² + v ² . Si s = 0 o s ≥ 1 , descartamos u y v , y probamos con otro par ( u , v ) . Debido a que u y v están distribuidos uniformemente y a que solo 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 comprobar 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${\textstyle {\sqrt {s}}}$, dividido por${\displaystyle \pi }$. A partir de esto, encontramos que la función de densidad de probabilidad tiene el valor constante 1 en el intervalo (0, 1) . Igualmente, el ángulo θ dividido por${\displaystyle 2\pi }$se distribuye uniformemente en el intervalo [ 0, 1) y es independiente de s .

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

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}}}}$$ y $${\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}}}.}$$

El método polar se diferencia del método básico en que es un tipo de muestreo por 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 costosas funciones trigonométricas mejora la velocidad con respecto a la forma básica. ^{[ 6 ]} Descarta 1 − π /4 ≈ 21,46% del total de pares de números aleatorios distribuidos uniformemente de entrada generados, es decir, descarta 4/ π − 1 ≈ 27,32% de pares de números aleatorios distribuidos uniformemente por cada par de números aleatorios gaussianos generados, lo que requiere 4/ π ≈ 1,2732 números aleatorios de entrada por cada 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 usando 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 desde C como una función intrínseca ) para calcular funciones complejas. $${\displaystyle \exp(iz)=e^{iz}=\cos z+i\sin z,\,}$$ y simplemente separa las partes reales de las imaginarias.

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

Dejar${\textstyle r={\sqrt {-\ln(u_{1})}}}$y${\textstyle z=2\pi u_{2}.}$Entonces $${\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.

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

El número positivo más pequeño que se puede representar en números de punto flotante de precisión simple y doble según el estándar IEEE es subnormal y es${\displaystyle 2^{-149}}$y${\displaystyle 2^{-1074}}$respectivamente. Los números normales más pequeños son${\displaystyle (2^{-126}-2^{-149})}$y${\displaystyle (2^{-1022}-2^{-1074})}$^{La mayoría de} los generadores de números aleatorios que producen un número de punto flotante con distribución uniforme no utilizan esta precisión, ya que funcionan simplemente convirtiendo un entero aleatorio de L bits a punto flotante y luego dividiéndolo por ^2L o 2L - 1. [ ^{11 ] Es} posible escribir un algoritmo que utilice esta precisión adicional, aunque para acceder a ella será necesario consumir más bits aleatorios. ^{[ 12 ]}

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

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