генерация пуассоновских переменных в С++

Из другой вопрос I спрашивали ранее, кажется, вы также можете приблизить poisson(750) к poisson(375) + poisson(375).

Если вы пойдете по пути "существующей библиотеки", ваш компилятор может уже поддерживать пакет C++11 std::random. Вот как вы его используете:

#include <random>
#include <ctime>
#include <iostream>

std::mt19937 mrand(std::time(0));  // seed however you want

typedef long unsigned int luint;

luint poisson(luint lambda)
{
    std::poisson_distribution<luint> d(lambda);
    return d(mrand);
}

int main()
{
    std::cout << poisson(750) << '\n';
    std::poisson_distribution<luint> d(750);
    std::cout << d(mrand) << '\n';
    std::cout << d(mrand) << '\n';
}

Я использовал его двумя способами выше:

1. Я пытался имитировать ваш существующий интерфейс.

2. Если вы создаете std::poisson_distribution со средним значением, более эффективно использовать это распределение снова и снова для одного и того же среднего значения (как это делается в main()).

Вот пример вывода для меня:

751
730
779

exp(-750) - очень маленькое число, очень близкое к наименьшему возможному двойному значению, поэтому ваша проблема числовая. В любом случае ваша сложность будет линейной по лямбда, поэтому алгоритм не очень эффективен для высоких лямбда. Если у вас нет веских причин кодировать это самостоятельно, использование существующей реализации библиотеки, вероятно, имеет смысл, поскольку эти числовые алгоритмы, как правило, чувствительны именно к проблемам точности, с которыми вы сталкиваетесь.

Поскольку вы используете только L в выражении (p>L), вы, по сути, тестируете (log(p) > -lambda). Это не очень полезное преобразование. Конечно, вам больше не нужен exp(-750), но вместо этого вы просто переполнитесь p.

Теперь p — это просто Π(mrand.rand()), а log(p) — это log(Π(mrand.rand())) есть Σ(log(mrand.rand()). Это дает вам необходимое преобразование:

double logp = 0;
do {
    k++;
    logp += log(mrand.rand());
} while( logp > -lambda);

double имеет только 11-битную экспоненту, но 52-битную мантиссу. Следовательно, это значительное увеличение численной стабильности. Плата за это заключается в том, что вам нужно log на каждой итерации вместо одной exp заранее.

В подобных ситуациях вам не нужно вызывать генератор случайных чисел более одного раза. Все, что вам нужно, это таблица кумулятивных вероятностей:

double c[k] = // the probability that X <= k (k = 0,...)

Затем сгенерируйте случайное число 0 <= r < 1 и возьмите первое целое число X такое, что c[X] > r. Вы можете найти этот X с помощью бинарного поиска.

Чтобы сгенерировать эту таблицу, нам нужны отдельные вероятности

p[k] = lambda^k / (k! e^lambda) // // the probability that X = k

Если lambda велико, это становится крайне неточным, как вы обнаружили. Но здесь мы можем использовать хитрость: начать с (или близкого) к наибольшему значению, с k = floor[lambda], и представить на данный момент, что p[k] равно 1. Затем вычислите p[i] для i > k, используя рекуррентное соотношение

p[i+1] = (p[i]*lambda) / (i+1)

и для i < k использования

p[i-1] = (p[i]*i)/lambda

Это гарантирует, что наибольшие вероятности имеют максимально возможную точность.

Теперь просто вычислите c[i], используя c[i+1] = c[i] + p[i+1], до точки, где c[i+1] совпадает с c[i]. Затем вы можете нормализовать массив, разделив на это предельное значение c[i]; или вы можете оставить массив как есть и использовать случайное число 0 <= r < c[i].

См.: http://en.wikipedia.org/wiki/Inverse_transform_sampling.