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

这几乎是 Knuth 的 TAOCP 第 3 版第 122 页的第 2 卷中算法的直接实现。

问题是rand_gauss()有时会返回区间 [0,1] 之外的值。

4

2 回答 2

8

Knuth 在 TAOCP 第 2 卷的第 122 页描述了极坐标法。该算法生成均值 = 0且标准差 = 1的正态分布。但是您可以通过乘以所需的标准偏差并添加所需的平均值来调整它。

您可能会发现将您的代码与 C-FAQ 中极坐标方法的另一个实现进行比较很有趣。

于 2011-03-13T03:40:02.070 回答
4

将您的 if 语句更改为(s >= 1.0 || s == 0.0). 更好的是，使用break以下示例中所示的 a 作为 SIMD 高斯随机数生成返回复数对 (u,v)。这使用了Mersenne twister 随机数生成器 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;
}

这个算法出奇的快。为四个不同的高斯随机数生成器计算两个随机数 (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() 的嵌套多项式实现将 Box-Muller 计算速度提高了 2.9 倍。SIMD Box Muller 和极坐标算法当然具有竞争力。它们也可以很容易地并行化。使用 gcc -Ofast -S，汇编代码转储显示平方根是 SIMD SSE2： sqrt --> sqrtsd %xmm0, %xmm0

评论：使用 gcc5 获得准确的时间真的很困难也很沮丧，但我认为这些都可以：截至 2016 年 2 月 3 日：DLW

[1] 相关链接：cython中的c malloc数组指针返回

[2] 算法比较，但不一定适用于 SIMD 版本： http: //www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

[3] 查尔斯·K·加勒特： http: //krisgarrett.net/papers/l2approx.pdf

于 2016-01-31T08:11:13.787 回答