在 C# 中实现 Box-Mueller 随机数生成器

Question

来自this question: Random number generator which gravitates numbers to any given number in range? 我做了一些研究，因为我以前遇到过这样的随机数生成器。我只记得“Mueller”这个名字，所以我想我在这里找到了它：

• Box-Mueller transform

我可以找到许多其他语言的实现，但我似乎无法在 C# 中正确实现。

这个页面，例如，The Box-Muller Method for Generating Gaussian Random Numbers 说代码应该是这样的（这不是 C#）：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

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 之间。很容易就变成1000，但是发生的概率极低。

Answer 2

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

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

Answer 3

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

Answer 4

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

另外，我认为这不是 Box-Mueller 变换，实际上是 Marsaglia 极坐标法。

Answer 5

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

既然它是统计的，为什么仅仅因为std.dev是1就硬限制在-1..1？任何一方都可以在统计上有所发挥，但仍能满足统计要求。