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

Я хотел реализовать простой метод выборки из полиномиального распределения на 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))

Я получаю много единиц (вероятность 62%), много двоек, немного троек и т. д.

Однако, когда я пробую приведенную ниже реализацию на C# (довольно простая выборка с обратным преобразованием для многочлена, основанная только на универсальной случайной величине), я получаю действительно странные результаты. Для всех 1000 образцов я часто нахожу все единицы. Иногда я нахожу все 3 (!!??). Результаты никогда не выглядят так, как вы ожидаете (и что вы получаете от функции python — попробуйте запустить ее самостоятельно несколько раз). Это действительно страшно, так как мы полагаемся на эти примитивы. Кто-нибудь знает, что может быть не так с версией С#?

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

random c# probability multinomial statistical-sampling
person Rohit Pandey    schedule 17.11.2017    source источник
comment
Вы инициализируете Random каждый раз, когда вызываете MultinomialSample. Если эти вызовы очень близки друг к другу, Random будет инициализирован с одним и тем же начальным числом (на основе системных часов). Попробуйте либо сделать Random приватным полем класса: private static Random random = new Random();, либо передать его в метод в качестве аргумента из Main, где оно будет инициализировано только один раз.   -  person Rufus L    schedule 18.11.2017

Ответы (2)


arrow_upward
2
arrow_downward

Вы инициализируете 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();
}
person Rufus L    schedule 17.11.2017

arrow_upward
1
arrow_downward

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

person Joe Blogs    schedule 12.11.2019
comment
Отличный улов. Спасибо! - person Rohit Pandey; 12.11.2019