这种在 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 设为私有 class 字段：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 的代码进行上述更改时，我得到了概率数组指定的所需频率分布。

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

What's wrong with this simple method to sample from multinomial in C#?

c#

random

probability

multinomial

statistical-sampling