Générateur de nombre aléatoire gaussien

Question

J'essaie de mettre en œuvre un générateur de nombres aléatoires distribués gaussien dans l'intervalle [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)); 
} 

Il est à peu près une mise en œuvre avant droite de l'algorithme en 2e volume de TAOCP 3ème édition page de Knuth 122.

Le problème est que rand_gauss() retourne parfois des valeurs en dehors de l'intervalle [0, 1].

Answer 1

Knuth décrit la méthode polaire sur p 122 du 2ème volume de TAOCP. Cet algorithme génère une distribution normale avec moyenne = 0 et écart-type = 1. Mais vous pouvez ajuster cela en multipliant par l'écart-type désiré et en ajoutant la moyenne désirée.

Vous pourriez trouver amusant de comparer votre code à une autre implémentation de the polar method in the C-FAQ.

Answer 2

Changez votre instruction if en (s >= 1.0 || s == 0.0). Mieux encore, utilisez un break comme indiqué dans l'exemple suivant pour un nombre aléatoire SIMD gaussien générant le renvoi d'une paire complexe (u, v). Ceci utilise le générateur de nombres aléatoires Mersennedsfmt(). Si vous ne voulez qu'un seul nombre aléatoire réel, renvoyez uniquement u et enregistrez le v pour la prochaine passe.

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

Cet algorithme est étonnamment rapide. Les temps d'exécution pour le calcul de deux nombres aléatoires (u, v) pour quatre générateurs de nombres aléatoires gaussiennes sont:

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) 

Les routines polynomiales fast_sin et fast_cos de Charles K. Garrett accélérer le calcul de Box-Muller d'un facteur 2,9 en utilisant une implémentation polynomiale imbriquée de cos() et sin(). Le SIMD Box Muller et les algorithmes polaires sont certainement compétitifs. En outre, ils peuvent être parallélisés facilement. En utilisant gcc -Ofast -S, le vidage du code d'assemblage montre que la racine carrée est la SIMD SSE2: sqrt -> sqrtsd% xmm0,% xmm0

Commentaire: il est vraiment difficile et frustrant d'obtenir des synchronisations précises avec gcc5, mais je pense que ce sont ok: au 03.02.2016: DLW

[1] lien connexe: c malloc array pointer return in cython

[2] Une comparaison des algorithmes, mais pas nécessairement pour les versions SIMD: http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

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