реализация распределения Пуассона в c ++

Я пытаюсь написать программу для вычисления вероятностной функции масс распределения Пуассона P (x = n) с параметром лямбда, используя эту формулу: ( (e^-lambda)*(lambda^n))/n!

Этот подход хорошо работает, когда я использую маленькие лямбды и маленькие числа, но если я хочу вычислить, например, P (x = 30) с лямбда 20, результат будет 4.68903e + 006, что неверно.

Думаю, проблема в вычислении n !. Я реализовал функцию для вычисления значения факториала и использовал тип данных unsigned long long для результата вычисления факториала, но проблема в том, что количество 30! равно 265,252,859,812,191,058,636,308,480,000,000, а максимальное число, доступное для длинных длин без знака, составляет 18,446,744,073,709,551,615, что меньше 30 !.

Что мне делать, чтобы решить эту проблему? Есть ли какой-либо другой способ или какая-либо функция для вычисления этой возможности в С ++?

тип данных


c++ distribution poisson factorial
person starrr    schedule 10.05.2015    source источник
comment
Почему голос против? Что не так с этим вопросом?   -  person Cheiron    schedule 11.05.2015
comment
@Cheiron: Трудно найти ошибки в коде, который мы не видим. Руководящие принципы SO довольно ясны по этому поводу.   -  person MSalters    schedule 11.05.2015
comment
@MSalters Вопросы не должны включать код, чтобы быть по теме SO. Фактически, многие вопросы с самым высоким рейтингом не имеют кода.   -  person isarandi    schedule 11.05.2015
comment
@isarandi: Я знаю, я слежу за Meta.SO. Но посмотрите на начало третьего абзаца: Думаю, проблема в .... Вот почему нам нужен код, если он существует.   -  person MSalters    schedule 11.05.2015
comment
Рекомендации SO также специфичны в отношении того факта, что вы должны комментировать, что не так при голосовании против сообщения.   -  person Cheiron    schedule 12.05.2015
comment
В статье Википедии о распределении Пуассона дается рекомендация реализовать его вероятностную массу функция (не путать с функцией случайной величины) в терминах функции std::lgamma, чтобы избежать проблем с большими числами.   -  person rwong    schedule 27.10.2015
comment
Если приближения в порядке, вы можете использовать распределение Гаусса, когда n равно ›15 или около того (или всякий раз, когда n! Переполняется).   -  person E L    schedule 13.12.2018

Ответы (3)


arrow_upward
3
arrow_downward

Одним из способов решения проблемы больших n является вычисление распределения в домене журнала:

X = ((e^-lambda)*(lambda^n))/n!
ln X = -lambda + n*ln(lambda) - Sum (ln(n))
return e^X
person Paulo Mendes    schedule 10.05.2015
comment
ITYM n * ln(lambda) - person MSalters; 11.05.2015
comment
ln (n!) = Сумма (ln (n))? - person alberto; 21.10.2017
comment
Да, где n внутри Sum фактически работает как индекс, а не как фиксированный n. Возможно, я слишком упростил его и, конечно, немного злоупотребил обозначениями. - person Paulo Mendes; 23.10.2017

arrow_upward
3
arrow_downward

Попробуй немного математики

 (lambda^n))/n!

Разве это не так

 (lambda/n) * (lambda/(n-1) * ...

и этими числами можно будет управлять с двойным, а не с очень большими числами.

person Ed Heal    schedule 10.05.2015
comment
Это плохая реализация. Оценка lambda^n - O (log N), оценка - O (N). Правильный способ реализовать это (в C ++ 98 до <random> дистрибутивов) - использовать приближение Стирлинга, который также является O (log N). (Еще лучший подход сначала вычисляет log(( (e^-lambda)*(lambda^n))/n!), поскольку это довольно тривиально) - person MSalters; 11.05.2015
comment
Это точно, а не приближение - person Ed Heal; 11.05.2015
comment
Это будет приближение в пределах трех сроков; компьютеры являются двоичными, поэтому lambda/n, lambda/n-1 или lambda/n-2 пытается разделить на три. - person MSalters; 11.05.2015
comment
Неужели O (n) при n = 30 действительно все так плохо в этом случае? Ответ прост. Единственное изменение, которое я мог бы сделать, это сделать следующее: double total = pow (lamda, n); for (int i = 2; i ‹= n; i ++) {total / = i; } - person E L; 13.12.2018

arrow_upward
3
arrow_downward

Если нужно только генерировать случайные значения из распределения Пуассона, а в противном случае не нужно знать его < em> функция массы вероятности, наиболее прямым способом является использование заранее заданного предопределенного дистрибутив, который является частью C ++ 11. Похожую реализацию также можно найти в Boost.

person MSalters    schedule 10.05.2015
comment
Я не знаю, как на самом деле извлечь значение pmf из std::poisson_distribution. Хотите уточнить? - person T.C.; 11.05.2015
comment
@ T.C. : Думаю, вы правы - не существует определенного сопоставления базового однородного диапазона с распределенным диапазоном Пуассона. Я собираюсь оставить здесь ответ на тот случай, если кто-то знает, как извлечь PMF из данного дистрибутива C ++. - person MSalters; 11.05.2015
comment
Я не думаю, что это возможно в общем случае. Во многих случаях вам не нужно вычислять PMF / PDF / CDF, чтобы сгенерировать число. Например, normal_distribution может использовать преобразование Бокса-Мюллера. - person T.C.; 13.05.2015