这种从 C# 中的多项式采样的简单方法有什么问题？

Question

我想在 C# 中实现一个从多项式分布中采样的简单方法（第一个参数是我们要采样的整数数组，第二个参数是选择每个整数的概率）。

当我在 python 中使用 numpy 执行此操作时，结果是有意义的。

np.random.choice(np.array([1,2,3,4,5,6]),p=np.array([.624,.23,.08,.04, .02, .006]),size=len(b))

我得到很多 1（概率 62%）、一堆 2、一些 3 等等。

但是，当我在 C# 中尝试下面的实现时（非常简单的多项式逆变换采样，仅依赖于一个统一的随机变量），我得到了非常奇怪的结果。对于所有 1000 个样本，我经常会找到全 1。有时，我会找到所有 3 的 (!!??)。结果永远不会像您期望的那样（以及您从 python 函数中获得的结果 - 尝试自己运行几次）。这真的很可怕，因为我们依赖这些原语。有没有人了解 C# 版本可能有什么问题？

    static void Main(string[] args)
    {
        int[] iis = new int[7];
        int[] itms = new int[] { 1, 2, 3, 4, 5, 6 };
        double[] probs = new double[] { .624, .23, .08, .04, .02, .006 };
        for (int i = 0; i < 1000; i++)
        {
            iis[MultinomialSample(itms, probs)] += 1;
        }

        foreach (var ii in iis)
        {
            Console.Write(ii + ",");
        }

        Console.Read();
    }


     private static int MultinomialSample(int[] s, double[] ps)
    {
        double[] cumProbs = new double[ps.Length];
        cumProbs[0] = ps[0];

        for (int i = 1; i < ps.Length; i++)
        {
            cumProbs[i] = cumProbs[i - 1] + ps[i];
        }

        Random random = new Random();
        double u = random.NextDouble();

        for (int i = 0; i < cumProbs.Length - 1; i++)
        {
            if (u < cumProbs[i])
            {
                return s[i];
            }
        }

        return s[s.Length - 1];
    }

Answer 1

每次调用 MultinomialSample 时，您都在初始化 Random。如果这些调用非常接近，Random 将使用相同的种子初始化（基于系统时钟）。尝试将 Random 设为私有类字段：private static Random random = new Random(); 或将其作为来自 Main 的参数传递给方法，在该方法中它只会被初始化一次：

private static Random random = new Random();

private static int MultinomialSample(IReadOnlyList<int> sample, 
    IReadOnlyList<double> probabilities)
{
    var cumProbs = new double[probabilities.Count];
    cumProbs[0] = probabilities[0];

    for (var i = 1; i < probabilities.Count; i++)
    {
        cumProbs[i] = cumProbs[i - 1] + probabilities[i];
    }

    for (var i = 0; i < cumProbs.Length - 1; i++)
    {
        if (random.NextDouble() < cumProbs[i])
        {
            return sample[i];
        }
    }

    return sample[sample.Count - 1];
}

private static void Main()
{
    var iis = new int[7];
    var items = new[] {1, 2, 3, 4, 5, 6};
    var probabilities = new[] {.624, .23, .08, .04, .02, .006};

    for (int i = 0; i < 1000; i++)
    {
        iis[MultinomialSample(items, probabilities)] ++;
    }

    Console.WriteLine(string.Join(", ", iis));

    Console.WriteLine("\nDone!\nPress any key to exit...");
    Console.ReadKey();
}

Answer 2

我在我正在研究的模拟中使用了 Rufus 的代码，并注意到仍然存在问题，即使在随机数生成器只播种一次之后（这是正确的做法）。您会注意到，在我们进行迭代时，对 random.NextDouble() 的调用每次都会生成一个新的随机数。这是错误的。

for (var i = 0; i < cumProbs.Length - 1; i++)
{
    if (random.NextDouble() < cumProbs[i])
    {
        return sample[i];
    }
}

随机数应该在循环外生成，如下：

var r = random.NextDouble();
for (var i = 0; i < cumProbs.Length - 1; i++)
{
    if (r < cumProbs[i])
    {
        return sample[i];
    }
}

您可以将其与维基百科上给出的 Excel 算法进行比较：https://en.wikipedia.org/wiki/Multinomial_distribution。当我对 Rufus 的代码进行上述更改时，我得到了概率数组指定的所需频率分布。