Что не так с этим простым методом для выборки из многочлена в C#?

Question

Я хотел реализовать простой метод для выборки из полиномиального распределения в C# (первый аргумент - это массив целых чисел, которые мы хотим выбрать, а второй - вероятности выбора каждого из этих целых чисел).

Когда я делаю это с Numpy в Python, результаты имеют смысл.

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

Вы инициализируете Random каждый раз, когда вы звоните MultinomialSample, Если эти звонки очень близко друг к другу, 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. Когда я внес вышеупомянутое изменение в код Руфуса, я получил желаемое частотное распределение, заданное массивом вероятностей.