用C＃实现Box-Mueller随机数发生器

从这个问题：随机数发生器将数字吸引到范围内的任何给定数字？ 我做过一些研究，因为我之前遇到过这样一个随机数发生器。 我记得的只是“穆勒”的名字，所以我想我找到了，在这里：

• Box-Mueller变换

我可以在其他语言中找到它的大量实现，但我似乎无法在C＃中正确实现它。

例如，这个页面，用于生成高斯随机数的Box-Muller方法表示代码应该如下所示（这不是C＃）：

#include  #include  #include  #include  double gaussian(void) { static double v, fac; static int phase = 0; double S, Z, U1, U2, u; if (phase) Z = v * fac; else { do { U1 = (double)rand() / RAND_MAX; U2 = (double)rand() / RAND_MAX; u = 2. * U1 - 1.; v = 2. * U2 - 1.; S = u * u + v * v; } while (S >= 1); fac = sqrt (-2. * log(S) / S); Z = u * fac; } phase = 1 - phase; return Z; } 

现在，这是我在C＃中实现的上述内容。 注意，变换产生2个数字，因此上面有“阶段”的技巧。 我只是丢弃第二个值并返回第一个值。

 public static double NextGaussianDouble(this Random r) { double u, v, S; do { u = 2.0 * r.NextDouble() - 1.0; v = 2.0 * r.NextDouble() - 1.0; S = u * u + v * v; } while (S >= 1.0); double fac = Math.Sqrt(-2.0 * Math.Log(S) / S); return u * fac; } 

我的问题是以下特定情况，我的代码没有返回0-1范围内的值，我无法理解原始代码如何。

• u = 0.5，v = 0.1
• S变为0.5*0.5 + 0.1*0.1 = 0.26
• fac变成~ 3.22
• 因此返回值为~ 0.5 * 3.22或~ 1.6

那不在0 .. 1 。

我做错了什么/不理解？

如果我修改我的代码，以便不是将fac乘以u ，而是乘以S ，我得到一个范围从0到1的值，但它的分布是错误的（似乎有一个大约0.7-0.8的最大分布，然后逐渐变细）在两个方向上。）

Answer 1

你的代码很好。 你的错误是认为它应该只在[0, 1]内返回值。 （标准）正态分布是在整个实线上具有非零权重的分布。 也就是说， [0, 1]之外的值是可能的。 实际上， [-1, 0]内的值与[-1, 0]内的值一样可能，而且[0, 1]的补码约为正态分布权重的66％。 因此，66％的时间我们期望值超出[0, 1] 。

此外，我认为这不是Box-Mueller变换，而实际上是Marsaglia极地方法。

Answer 2

我不是数学家，也不是统计学家，但如果我想到这一点，我不会指望高斯分布在精确范围内返回数字。 给定您的实现，平均值为0，标准偏差为1，因此我希望值在钟形曲线上分布，中心为0，然后随着数字偏离任何一侧的0而减小。 所以序列肯定会涵盖两个+/-数字。

那么因为它是统计的，为什么仅仅因为std.dev是1而难以限制在-1..1？ 在统计上可以在任何一方发挥一些作用，仍然满足统计要求。

Answer 3

均匀随机变量确实在0..1之内，但高斯随机变量（这是Box-Muller算法生成的）可以是实线上的任何位置。 有关详细信息，请参阅wiki / NormalDistribution 。

Answer 4

我认为该函数返回极坐标。 因此，您需要两个值才能获得正确的结果。

此外，高斯分布不在0 .. 1之间。 它很容易以1000结尾，但这种发生的可能性极低。

Answer 5

这是一个蒙特卡罗方法，因此您无法限制结果，但您可以做的是忽略样本。

 // return random value in the range [0,1]. double gaussian_random() { double sigma = 1.0/8.0; // or whatever works. while ( 1 ) { double z = gaussian() * sigma + 0.5; if (z >= 0.0 && z <= 1.0) return z; } }