генерация пуассоновых переменных в с ++ - PullRequest
Новая
       103

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

5 голосов
/ 14 апреля 2011

Я реализовал эту функцию для генерации пуассоновской случайной величины 

typedef long unsigned int luint;
luint poisson(luint lambda) {
    double L = exp(-double(lambda));
    luint k = 0;
    double p = 1;
    do {
        k++;
        p *= mrand.rand();
    } while( p > L);
    return (k-1);
}

где mrand - генератор случайных чисел MersenneTwister. Я считаю, что при увеличении лямбды ожидаемое распределение будет неправильным со средним значением, которое насыщается на уровне около 750. Это из-за числовых приближений или я допустил какие-либо ошибки?

Ответы [ 5 ]

2 голосов
/ 14 апреля 2011

Если вы идете по пути "существующей библиотеки", ваш компилятор может уже поддерживать пакет 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
2 голосов
/ 14 апреля 2011

Поскольку вы используете только L в выражении (p>L), вы проверяете (log(p) > -lambda). Это не очень полезное преобразование. Конечно, вам больше не нужен опыт (-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 заранее.

2 голосов
/ 14 апреля 2011

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

1 голос
/ 14 апреля 2011

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

0 голосов
/ 14 апреля 2011

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

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

...