Método polar de Marsaglia

Método para generar números pseudoaleatorios

El método polar de Marsaglia ^[1] es un método de muestreo de números pseudoaleatorios para generar un par de variables aleatorias normales estándar independientes . ^[2]

Las variables aleatorias normales estándar se utilizan con frecuencia en informática , estadística computacional y, en particular, en aplicaciones del método de Monte Carlo .

El método polar funciona eligiendo puntos aleatorios ( x ,  y ) en el cuadrado −1 <  x  < 1, −1 <  y  < 1 hasta que

${\displaystyle 0<s=x^{2}+y^{2}<1,\,}$

y luego devolver el par requerido de variables aleatorias normales como

${\displaystyle x{\sqrt {\frac {-2\ln(s)}{s}}}\,,\ \ y{\sqrt {\frac {-2\ln(s)}{s}}},}$

o, equivalentemente,

${\displaystyle {\frac {x}{\sqrt {s}}}{\sqrt {-2\ln(s)}}\,,\ \ {\frac {y}{\sqrt {s}}}{\sqrt {-2\ln(s)}},}$

donde y representan el coseno y el seno del ángulo que forma el vector ( x , y ) con el eje x . ${\displaystyle x/{\sqrt {s}}}$ ${\displaystyle y/{\sqrt {s}}}$

Fundamento teórico

La teoría subyacente puede resumirse de la siguiente manera:

Si u se distribuye uniformemente en el intervalo 0 ≤  u  < 1, entonces el punto (cos(2π u ), sin(2π u )) se distribuye uniformemente en la circunferencia unitaria x ²  +  y ²  = 1, y multiplicando ese punto por una variable aleatoria independiente ρ cuya distribución es

${\displaystyle \Pr(\rho <a)=\int _{0}^{a}re^{-r^{2}/2}\,dr}$

producirá un punto

${\displaystyle \left(\rho \cos(2\pi u),\rho \sin(2\pi u)\right)}$

cuyas coordenadas se distribuyen conjuntamente como dos variables aleatorias normales estándar independientes.

Historia

Esta idea se remonta a Laplace , a quien Gauss atribuye el hallazgo de la ecuación anterior.

${\displaystyle I=\int _{-\infty }^{\infty }e^{-x^{2}/2}\,dx}$

tomando la raíz cuadrada de

${\displaystyle I^{2}=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }e^{-(x^{2}+y^{2})/2}\,dx\,dy=\int _{0}^{2\pi }\int _{0}^{\infty }re^{-r^{2}/2}\,dr\,d\theta .}$

La transformación a coordenadas polares evidencia que θ está uniformemente distribuida (densidad constante) de 0 a 2π, y que la distancia radial r tiene densidad

${\displaystyle re^{-r^{2}/2}.\,}$

( r ^{2 tiene la distribución} de chi cuadrado apropiada ).

Este método de producir un par de variables normales estándar independientes mediante la proyección radial de un punto aleatorio en la circunferencia unitaria a una distancia dada por la raíz cuadrada de una variable chi-cuadrado-2 se denomina método polar para generar un par de variables aleatorias normales.

Consideraciones prácticas

Una aplicación directa de esta idea,

${\displaystyle x={\sqrt {-2\ln(u_{1})}}\cos(2\pi u_{2}),\quad y={\sqrt {-2\ln(u_{1})}}\sin(2\pi u_{2})}$

se llama transformada de Box-Muller , en la que la variable chi se genera generalmente como

${\displaystyle {\sqrt {-2\ln(u_{1})}},}$

pero esa transformación requiere funciones de logaritmo, raíz cuadrada, seno y coseno. En algunos procesadores, el coseno y el seno del mismo argumento se pueden calcular en paralelo utilizando una única instrucción. ^[3] En particular, para las máquinas basadas en Intel, se puede utilizar la instrucción de ensamblador fsincos o la instrucción expi (disponible, por ejemplo, en D) para calcular ecuaciones complejas.

${\displaystyle \operatorname {expi} (z)=e^{iz}=\cos(z)+i\sin(z),\,}$

y simplemente separar las partes reales e imaginarias.

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

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

${\displaystyle \ re^{iz}={\sqrt {-2\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].}$

En cambio, el método polar elimina la necesidad de calcular un coseno y un seno. En cambio, al resolver un punto en el círculo unitario, estas dos funciones se pueden reemplazar con las coordenadas x e y normalizadas al radio. En particular, un punto aleatorio ( x ,  y ) dentro del círculo unitario se proyecta sobre la circunferencia unitaria al establecer y formar el punto ${\displaystyle {\sqrt {x^{2}+y^{2}}}}$ ${\displaystyle s=x^{2}+y^{2}}$

${\displaystyle \left({\frac {x}{\sqrt {s}}},{\frac {y}{\sqrt {s}}}\right),\,}$

que es un procedimiento más rápido que calcular el coseno y el seno. Algunos investigadores sostienen que la instrucción condicional if (para rechazar un punto fuera del círculo unitario) puede hacer que los programas sean más lentos en los procesadores modernos equipados con segmentación y predicción de saltos. ^[4] Además, este procedimiento requiere aproximadamente un 27 % más de evaluaciones del generador de números aleatorios subyacente (solo los puntos generados se encuentran dentro del círculo unitario). ${\displaystyle \pi /4\approx 79\%}$

Ese punto aleatorio en la circunferencia se proyecta entonces radialmente a la distancia aleatoria requerida por medio de

${\displaystyle {\sqrt {-2\ln(s)}},\,}$

utilizando la misma s porque s es independiente del punto aleatorio en la circunferencia y está distribuida uniformemente de 0 a 1.

Implementación

Java

Implementación simple en Java utilizando la media y la desviación estándar:

privado estático doble repuesto ; privado estático booleano hasSpare = falso ;        public static sincronizado doble generateGaussian ( doble media , doble stdDev ) { if ( hasSpare ) { hasSpare = false ; return spare * stdDev + media ; } else { double u , v , s ; do { u = Math.random ( ) * 2-1 ; v = Math.random ( ) * 2-1 ; s = u * u + v * v ; } while ( s > = 1 || s == 0 ) ; s = Math.sqrt ( -2.0 * Math.log ( s ) / s ) ; spare = v * s ; hasSpare = true ; return mean + stdDev * u * s ; } }​​​​​​​​                                                                                     

C++

Una implementación no segura para subprocesos en C++ que utiliza la media y la desviación estándar:

doble generateGaussian ( doble media , doble desviación estándar ) { static double spare ; static bool hasSpare = false ;              si ( tiene repuesto ) { tiene repuesto = falso ; devolver repuesto * desviación estándar + media ; } de lo contrario { doble u , v , s ; hacer { u = ( rand () / (( doble ) RAND_MAX )) * 2,0 - 1,0 ; v = ( rand () / (( doble ) RAND_MAX )) * 2,0 - 1,0 ; s = u * u + v * v ; } mientras ( s >= 1,0 || s == 0,0 ); s = sqrt ( -2,0 * log ( s ) / s ); repuesto = v * s ; tiene repuesto = verdadero ; devolver media + desviación estándar * u * s ; } }                                                                                

La implementación de std::normal_distribution en C++11 GNU GCC libstdc++ utiliza el método polar de Marsaglia, como se cita aquí.

Julia

Una implementación sencilla de Julia :

"""  muestra de marsaglia(N)Genere `2N` muestras de la distribución normal estándar utilizando el método de Marsaglia. """ function marsagliasample ( N ) z = Array { Float64 }( undef , N , 2 ); for i in axes ( z , 1 ) s = Inf ; while s > 1 z [ i , : ] . = 2 rand ( 2 ) .- 1 ; s = sum ( abs2 . ( z [ i ,: ])) end z [ i , : ] .*= sqrt ( -2 log ( s ) / s ) ; end vec ( z ) end                             """  muestra de marsaglia(n,μ,σ)Generar `n` muestras de la distribución normal con media `μ` y desviación estándar `σ` utilizando el método de Marsaglia. """ function marsagliasample ( n , μ , σ ) μ .+ σ * marsagliasample ( cld ( n , 2 ))[ 1 : n ]; fin    

El bucle for se puede paralelizar mediante el uso de la Threads.@threadsmacro.

Referencias

1. Marsaglia, G.; Bray, TA (1964). "Un método conveniente para generar variables normales". SIAM Review . 6 (3): 260–264. Bibcode :1964SIAMR...6..260M. doi :10.1137/1006063. JSTOR  2027592.
2. Peter E. Kloeden Eckhard Platen Henri Schurz, Solución numérica de SDE mediante experimentos informáticos, Springer, 1994.
3. Kanter, David. "Arquitectura gráfica Ivy Bridge de Intel". Real World Tech . Consultado el 8 de abril de 2013 .
4. Este efecto puede verse acentuado en una GPU que genera muchas variables en paralelo, donde un rechazo en un procesador puede ralentizar a muchos otros procesadores. Véase la sección 7 de Thomas, David B.; Howes, Lee W.; Luk, Wayne (2009), "A comparison of CPUs, GPUs, FPGAs, and massively parallel processing arrays for random number generation", en Chow, Paul; Cheung, Peter YK (eds.), Proceedings of the ACM/SIGDA 17th International Symposium on Field Programmable Gate Arrays, FPGA 2009, Monterey, California, EE. UU., 22-24 de febrero de 2009 , Association for Computing Machinery , pp. 63-72, CiteSeerX 10.1.1.149.6066 , doi :10.1145/1508128.1508139, ISBN  9781605584102, Número de identificación del sujeto  465785.