8

区間[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));
}

これは、クヌースのTAOCP第3版122ページの第2巻にあるアルゴリズムの非常に単純な実装です。

問題は、rand_gauss（）が間隔[0,1]外の値を返すことがあることです。

4

2 に答える 2

8

クヌースは、TAOCPの第2巻のp122で極座標法について説明しています。このアルゴリズムは、平均=0および標準偏差=1の正規分布を生成します。ただし、必要な標準偏差を掛けて、必要な平均を加算することにより、これを調整できます。

コードをC-FAQのpolarメソッドの別の実装と比較するのは楽しいかもしれません。

于 2011-03-13T03:40:02.070 に答える
4

ifステートメントをに変更します(s >= 1.0 || s == 0.0)。さらに良いbreakことに、次の例に示すように、複素数のペア（u、v）を返すSIMDガウス乱数を生成するためにを使用します。これは、メルセンヌツイスター乱数ジェネレーター dsfmt()を使用します。単一の実数の乱数のみが必要な場合は、戻り値のみを返し、次のパスuのためにを保存します。v

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

このアルゴリズムは驚くほど高速です。4つの異なるガウス乱数ジェネレーターの2つの乱数（u、v）を計算するための実行時間は次のとおりです。

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)

Charles K. Garrettのfast_sinおよびfast_cos多項式ルーチンは、cos（）およびsin（）のネストされた多項式実装を使用して、ボックスミュラー計算を2.9倍高速化します。SIMDボックスミュラー法と極座標アルゴリズムは確かに競争力があります。また、それらは簡単に並列化できます。gcc -Ofast -Sを使用すると、アセンブリコードダンプは平方根がSIMD SSE2であることを示します：sqrt-> sqrtsd％xmm0、％xmm0

コメント：gcc5で正確なタイミングを取得するのは本当に難しくてイライラしますが、これらは問題ないと思います：2016年2月3日現在：DLW

[1]関連リンク：cmalloc配列ポインターがcythonで返される

[2]アルゴリズムの比較。ただし、SIMDバージョンである必要はありません：http：//www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

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

于 2016-01-31T08:11:13.787 に答える