Generador de números aleatorios de Gauss

Question

Estoy tratando de implementar un generador de números aleatorios distribuidos de Gauss en el intervalo [0,1].

float rand_gauss (void) { 
    float v1,v2,s; 

    do { 
    v1 = 2.0 * ((float) rand()/RAND_MAX) - 1; 
    v2 = 2.0 * ((float) rand()/RAND_MAX) - 1; 

    s = v1*v1 + v2*v2; 
    } while (s >= 1.0); 

    if (s == 0.0) 
    return 0.0; 
    else 
    return (v1*sqrt(-2.0 * log(s)/s)); 
} 

Es más o menos una aplicación recta hacia adelante del algoritmo en segundo volumen de TAOCP tercera página de edición de Knuth 122.

El problema es que rand_gauss() veces devuelve valores fuera del intervalo [0, 1].

Answer 1

Knuth describe el método polar en p 122 del 2 ° volumen de TAOCP. Ese algoritmo genera una distribución normal con mean = 0 y desviación estándar = 1. Pero puede ajustar eso multiplicando por la desviación estándar deseada y agregando la media deseada.

Puede que le resulte divertido comparar su código con otra implementación de the polar method in the C-FAQ.

Answer 2

Cambie su estado de cuenta if al (s >= 1.0 || s == 0.0). Mejor aún, use un break como se ve en el siguiente ejemplo para un número aleatorio gaussiano SIMD que genera la devolución de un par complejo (u, v). Esto utiliza el generador de números aleatorios twister Mersennedsfmt(). Si solo desea un número aleatorio único, real, devuelva solo u y guarde el v para la siguiente pasada.

inline static void randn(double *u, double *v) 
{ 
double s, x, y; // SIMD Marsaglia polar version for complex u and v 
while (1){ 
    x = dsfmt_genrand_close_open(&dsfmt) - 1.; 
    y = dsfmt_genrand_close_open(&dsfmt) - 1.; 
    s = x*x + y*y; 
    if (s < 1) break; 
} 
s = sqrt(-2.0*log(s)/s); 
*u = x*s; *v = y*s; 
return; 
} 

Este algoritmo es sorprendentemente rápido. Los tiempos de ejecución para el cálculo de dos números al azar (u, v) durante cuatro generadores de números aleatorios gaussianos diferentes son:

Times for delivering two Gaussian numbers (u + iv) 
i7-2600K @ 4GHz, gcc -Wall -Ofast -msse2 .. 

gsl_ziggurat      =  20.3 (ns) 
Box-Muller       =  78.8 (ns) 
Box-Muller with fast_sin fast_cos =  28.1 (ns) 
SIMD Marsaglia polar    =  35.0 (ns) 

El fast_sin y fast_cos rutinas polinómicas de Charles K. Garrett acelerar el cálculo de Box-Muller por un factor de 2,9 utilizando una implementación polinómica anidada de cos() y sin(). SIMD Box Muller y los algoritmos polares son ciertamente competitivos. También se pueden paralelizar fácilmente. Usando gcc -Ofast -S, el volcado de código ensamblador muestra que la raíz cuadrada es SIMD SSE2: sqrt -> sqrtsd% xmm0,% xmm0

Comentario: es realmente difícil y frustrante obtener sincronizaciones precisas con gcc5, pero creo que estos son bien: a partir del 02/03/2016: DLW

[1] enlace relacionado: c malloc array pointer return in cython

[2] Una comparación de algoritmos, pero no necesariamente para las versiones SIMD: http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

[3] Charles K. Garrett: http://krisgarrett.net/papers/l2approx.pdf