Transformada de Box-Muller

Visualización de la caja – Muller transforman – los puntos coloreados en la plaza unidad (u1, u2), dibujados como círculos, se mapean a un Gausiano 2D (z0, z1), dibujado como cruces. Las parcelas en los márgenes son las funciones de distribución de probabilidad de z0 y z1. Tenga en cuenta que 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, pasar por encima de un punto para destacarlo y su punto correspondiente.

La transformada Box-Muller, de George Edward Pelham Box y Mervin Edgar Muller, 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) números, dada una fuente de números aleatorios uniformemente distribuidos. De hecho, el método fue mencionado explícitamente por primera vez por Raymond E. A. C. Paley y Norbert Wiener en 1934.

La transformada de Box-Muller se suele expresar 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 de 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. El algoritmo de ziggurat brinda 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).

Forma básica

Suponga que U₁ y U₂ son muestras independientes elegidas de la distribución uniforme en el intervalo unitario (0, 1). Dejar

Z0=R#⁡ ⁡ ().. )=− − 2In⁡ ⁡ U1#⁡ ⁡ ()2π π U2){displaystyle ¿Qué?

y

Z1=Rpecado⁡ ⁡ ().. )=− − 2In⁡ ⁡ U1pecado⁡ ⁡ ()2π π U2).{displaystyle Z_{1}=Rsin(Theta)={sqrt {-2ln U_{1}sin(2pi U_{2}).,}

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

La derivación se basa en una propiedad de un sistema cartesiano bidimensional, donde las coordenadas X e Y están descritas por dos variables aleatorias independientes y normalmente distribuidas, las variables aleatorias para R² y Θ (que se muestra arriba) en las coordenadas polares correspondientes también son independientes y se pueden expresar como

R2=− − 2⋅ ⋅ In⁡ ⁡ U1{displaystyle R^{2}=-2cdot ln U_{1},}

y

.. =2π π U2.{displaystyle Theta =2pi U_{2}.

Porque R² 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 de chi-cuadrado coincide con la distribución exponencial, y la ecuación para R² anterior es una forma simple de generar el requerido variable exponencial.

Forma polar

Dos valores distribuidos uniformemente, u y v se utilizan para producir el valor s = R², que también se distribuye uniformemente. Las definiciones del seno y el cosino se aplican luego a la forma básica de la transformación Box-Muller para evitar el uso de funciones trigonométricas.

La forma polar fue propuesta por primera vez por J. Bell y luego modificada por R. Knop. 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 Recetas numéricas.

Dado u y v, independiente y distribuido uniformemente en el intervalo cerrado [−1, +1], conjunto s = R² = u² + v². Si s = 0 o s ≥ 1, descarte u y v, y probar otro par (u,v). Porque... u y v se distribuyen uniformemente y porque sólo se han admitido puntos dentro del círculo unitario, los valores de s se distribuirá uniformemente en el intervalo abierto (0, 1), también. Este último se puede ver calculando la función de distribución acumulativa para s en el intervalo (0, 1). Esta es la zona de un círculo con radio s{displaystyle scriptstyle {sqrt {s}}, dividida por π π {displaystyle scriptstyle pi }. De esto encontramos la función de densidad de probabilidad para tener el valor constante 1 en el intervalo (0, 1). Igualmente, el ángulo θ dividido por 2π π {displaystyle scriptstyle 2pi} se distribuye uniformemente en el intervalo [0, 1) e independiente de s.

Ahora identificamos el valor de s con la de U₁ y Silencio Silencio /()2π π ){displaystyle scriptstyle theta /(2pi)} con la de U₂ en la forma básica. Como se muestra en la figura, los valores de #⁡ ⁡ Silencio Silencio =#⁡ ⁡ 2π π U2{displaystyle scriptstyle cos theta =cos 2pi U_{2} y pecado⁡ ⁡ Silencio Silencio =pecado⁡ ⁡ 2π π U2{displaystyle scriptstyle sin theta =sin 2pi U_{2} en la forma básica se puede reemplazar con las ratios #⁡ ⁡ Silencio Silencio =u/R=u/s{displaystyle scriptstyle cos theta =u/R=u/{sqrt {s}} y pecado⁡ ⁡ Silencio Silencio =v/R=v/s{displaystyle scriptstyle sin theta =v/R=v/{sqrt {s}}, respectivamente. La ventaja es que calcular las funciones trigonométricas directamente se puede evitar. Esto es útil cuando las funciones trigonométricas son más costosas de calcular que la única división que reemplaza a cada una.

Así como la forma básica produce dos desviaciones normales estándar, también lo hace este cálculo alternativo.

z0=− − 2In⁡ ⁡ U1#⁡ ⁡ ()2π π U2)=− − 2In⁡ ⁡ s()us)=u⋅ ⋅ − − 2In⁡ ⁡ ss{displaystyle {fn0}csqrt {-2ln U_{1}}cos(2pi U_{2}={sqrt {-2ln s}left({frac {u}sqrt {s}right)=ucdot {fnMicroc} {-2ln s} {}}}}

y

z1=− − 2In⁡ ⁡ U1pecado⁡ ⁡ ()2π π U2)=− − 2In⁡ ⁡ s()vs)=v⋅ ⋅ − − 2In⁡ ⁡ ss.{displaystyle z_{1}={sqrt {-2ln U_{1}sin(2pi U_{2}={sqrt {-2ln s}left({frac {sqrt {}right)=vcdot {fnMicroc} {-2ln s} {}}}

Contraste de las dos formas

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 simple de calcular (siempre que el generador de números aleatorios sea relativamente rápido) y es numéricamente más robusto. Evitar el uso de costosas funciones trigonométricas mejora la velocidad con respecto a la forma básica. Descarta 1 − π/4 ≈ 21,46% de la entrada total de manera uniforme pares de números aleatorios distribuidos generados, es decir, descarta 4/π − 1 ≈ 27,32 % pares de números aleatorios distribuidos uniformemente por par de números aleatorios gaussianos generados, 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. En algunos procesadores, el coseno y el seno del mismo argumento se pueden calcular en paralelo usando una sola instrucción. En particular, para las máquinas basadas en Intel, se puede usar la instrucción del ensamblador fsincos o la instrucción expi (generalmente disponible en C como una función intrínseca), para calcular complejos

exp⁡ ⁡ ()iz)=eiz=#⁡ ⁡ z+ipecado⁡ ⁡ z,{displaystyle exp(iz)=e^{iz}=cos z+isin z,

y simplemente separa las partes real e imaginaria.

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

Vamos r=− − 2In⁡ ⁡ ()u1){displaystyle r={sqrt {-2ln(u_{1}}}} y z=2π π u2.{displaystyle z=2pi u_{2} Entonces...

reiz=− − 2In⁡ ⁡ ()u1)ei2π π u2=− − 2In⁡ ⁡ ()u1)[#⁡ ⁡ ()2π π u2)+ipecado⁡ ⁡ ()2π π u2)].{displaystyle {fnK}}e^{i2pi) u_{2}={sqrt {-2ln(u_{1}}}left[cos(2pi u_{2})+isin(2pi 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 una computadora se utiliza para producir una variable aleatoria uniforme tendrá inevitablemente algunas imprecisiones porque hay un límite inferior en cuanto los números cercanos pueden ser a 0. Si el generador utiliza 32 bits por valor de salida, el menor número no cero que se puede generar es 2− − 32{displaystyle 2^{-32}. Cuando U1{displaystyle U_{1} y U2{displaystyle U_{2} son iguales a esto que la transformación Box-Muller produce una desviación normal al azar igual a δ δ =− − 2In⁡ ⁡ ()2− − 32)#⁡ ⁡ ()2π π 2− − 32).. 6.660{displaystyle delta ={sqrt {-2ln(2^{-32}}cos(2pi 2^{-32})approx 6.660}. 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 2()1− − CCPR CCPR ()δ δ ))≃ ≃ 2.738× × 10− − 11{displaystyle 2(1-Phi (delta))simeq 2.738times 10^{-11} perdido por la truncación, donde CCPR CCPR ()δ δ ){displaystyle Phi (delta)} es la distribución normal acumulativa estándar. Con 64 bits el límite es empujado a δ δ =9.419{displaystyle delta =9.419} desviaciones estándar, para las cuales <math alttext="{displaystyle 2(1-Phi (delta))2()1− − CCPR CCPR ()δ δ )).5× × 10− − 21{displaystyle 2(1-Phi (delta)))Seleccionado5times 10^{-21}<img alt="{displaystyle 2(1-Phi (delta)).

Implementación

El cuadro estándar – Muller transforma genera valores de la distribución normal estándar (i.e. desviaciones normales) con media 0 y desviación estándar 1. La implementación de abajo en la norma C+ genera valores de cualquier distribución normal con media μ μ {displaystyle mu } y diferencia σ σ 2{displaystyle sigma ^{2}. Si Z{displaystyle Z} es una desviación normal normal, entonces X=Zσ σ +μ μ {displaystyle X=Zsigma +mu} tendrá una distribución normal con media μ μ {displaystyle mu } y desviación estándar σ σ {displaystyle sigma }. Tenga en cuenta que el generador de números aleatorios ha sido sembrado para asegurar que los nuevos valores de pseudo-aleatoria serán devueltos de llamadas secuenciales a los generateGaussianNoise función.

#include - No.#include > >#include - No.#include Identificadastd::par.doble, doble■ generarGaussianNoise()doble mu, doble sigma){} constexpr doble epsilon = std::numeric_limits.doble- No.epsilon(); constexpr doble 2_pi = 2.0 * M_PI; //inicializar el generador de números uniformes al azar (runif) en un rango 0 a 1 estática std::mt19937 rng()std::random_device{}()); // Standard mersenne_twister_engine seeded with rd() estática std::uniform_real_distribución■ Runif()0,0, 1.0); //crear dos números aleatorios, asegúrese de que u1 es mayor que epsilon doble u1, u2; do {} u1 = Runif()rng); } mientras ()u1 . epsilon); u2 = Runif()rng); //compute z0 y z1 auto mag = sigma * Sqrt()-2.0 * log()u1)); auto z0 = mag * #()2_pi * u2) + mu; auto z1 = mag * pecado()2_pi * u2) + mu; retorno std::make_pair()z0, z1);}