Gaussien générateur de nombre aléatoire

Je suis en train de mettre en œuvre une gaussienne distribué générateur de nombre aléatoire 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));
}

C'est à peu près tout droit vers l'avant de la mise en œuvre de l'algorithme de Knuth du 2e volume de TAOCP 3e édition, page 122.

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

Une gaussienne est sans limite. Ai-je raté quelque chose?
il y a une variance et de la moyenne, je prends la moyenne 0 et de variance^2 comme 1, de la distribution normale standard.
Une distribution normale standard peut prendre n'importe quelle valeur entre -l'infini et l'infini, avec une certaine probabilité; il n'y a pas de limite sur le résultat.
Est-il une raison pour laquelle vous êtes à l'aide de float plutôt que double? Habituellement, c'est une mauvaise idée..
L'arithmétique sur float et double est presque certainement le même coût, de plus, vous êtes à la conversion d'avant en arrière pour double de toute façon lorsque vous appelez log et sqrt.

OriginalL'auteur | 2011-03-13

  1. 8

    Knuth décrit le polar méthode à la page 122 du 2e volume de TAOCP. Cet algorithme génère une distribution normale avec moyenne = 0 et écart-type = 1. Mais vous pouvez régler que par la multiplication par souhaitée de l'écart type et l'ajout souhaité en dire.

    Vous trouverez peut-être amusant de comparer votre code à l'autre de la mise en œuvre de le polar à la méthode C-FAQ.

    OriginalL'auteur Mike Sherrill 'Cat Recall'

  2. 4

    Changer votre instruction if pour (s >= 1.0 || s == 0.0). Mieux encore, utilisez un break comme on le voit dans l'exemple suivant pour une SIMD Gaussien nombre aléatoire générant retour d'un complexe paire (u,v). Il utilise la Mersenne twister générateur de nombre aléatoire dsfmt(). Si vous ne voulez qu'une seule, la vraie, de nombres aléatoires, de retour seulement u et enregistrer la v pour le prochain passage.

    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 les quatre différentes Gaussien générateurs de nombres aléatoires 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)

    La fast_sin et fast_cos polynôme routines de Charles K. Garrett accélérer le Box-Muller calcul par un facteur de 2,9 à l'aide d'un imbriquée polynôme de mise en œuvre de cos() et sin(). Le SIMD Boîte de Muller et polaires algorithmes sont certainement concurrentiel. Ils peuvent aussi être facilement parallélisable. En utilisant gcc -Ofast -S, l'assemblée de code dump montre que la racine carrée est le SIMD SSE2: sqrt --> sqrtsd %xmm0, %xmm0

    Commentaire: il est vraiment difficile et frustrant d'obtenir exacts timings avec gcc5, mais je pense que ce sont ok: comme de 2/3/2016: DLW

    [1] lien: c malloc pointeur sur le tableau de retour en cython

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

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

    yup, merci, a changé. 73
    ici est la partie la plus difficile d'obtenir le droit: gcc -Wall-Ofast -msse2 -frename-registres -malign-double -fno-strict-aliasing DHAVE_SSE2=1 -DDSFMT_MEXP=19937 -o randn_test dSFMT.c randn_test.c -lm -
    Je suppose que @W9DKI, vous avez créé un double compte par erreur. Lorsque vous vous connectez la prochaine fois, de connexion de la même manière que vous avez fait la première façon et pas autrement. Aussi le drapeau de la fusion de comptes. Voir ce stackoverflow.com/help/merging-accounts
    Jonathan, merci pour la modification de la s=0 typo pour s==0. Au lieu d'utiliser while(1), un for(;;) serait un peu plus serré, mais gcc-Ofast compile à la fois pour la même chose.
    il est de toute évidence un facteur de 2 à gauche, c'est à dire x = 2.*dsfmt_genrand_close_open(&dsfmt) - 1.0; y = 2.*dsfmt_genrand_close_open(&dsfmt) - 1.0

    OriginalL'auteur W9DKI