Простое приближение кумулятивного распределения Пуассона в длинном хвосте?

10

Я хочу определить емкость C таблицы так, чтобы она имела остаточные шансы менее 2p переполнить для данного p[40120] , предполагая, что число записей следует закону Пуассона с заданной ожидаемой величиной .E[1031012]

В идеале я хочу, чтобы наименьшее целое число было Cтаким, что 1-CDF[PoissonDistribution[E],C] < 2^-pдля заданного pи E; но я доволен Cнемного выше этого. Mathematica отлично подходит для ручных вычислений, но я хотел бы вычислять Cс pи Eво время компиляции, что ограничивает меня 64-битной целочисленной арифметикой.

Обновление: в Mathematica (версия 7) e = 1000; p = 40; c = Quantile[PoissonDistribution[e], 1 - 2^-p]все 1231выглядит правильно (спасибо @Procrastinator); однако результат для обоих p = 50и p = 60есть 1250, что неправильно с небезопасной стороны (и имеет значение: мой эксперимент повторяется примерно в раз или более, и я хочу явно меньше, чем общих шансов отказа). Мне нужна грубая, но безопасная аппроксимация с использованием только 64-битной целочисленной арифметики , доступной в C (++) во время компиляции. 2 - 30225230

fgrieu
источник
1
Как насчет C = Quantile[PoissonDistribution[E],1-2^p]?
1
Главный член функции вероятности массы Пуассона доминирует в хвосте.
кардинал
1
@Procrastinator: да, это работает в Mathematica (за исключением признаков p, проблем точности, имен Eи Cзарезервированных). НО мне нужно простое приближение к этому, возможно, грубое (но безопасное) использование только 64-битной целочисленной арифметики!
fgrieu
3
По поводу обновления: Mathematica 8 возвращает 1262 для и 1290 для p = 60 . Нормальное приближение (@Proc): нельзя ожидать, что это будет хорошо работать в хвостах, что имеет решающее значение для расчета. p=50p=60
whuber
1
Возможно, вам следует спросить о stackoverflow. Я не знаком с вашими ограничениями. Я не знаю, что мешает вам использовать динамическое выделение памяти, или вы можете использовать ветвление для определения размера массива, или какова стоимость определения массива, который в два раза больше необходимого вам размера (а затем не использует все этого). Если какая-то функция, такая как (просто в качестве примера) дал вам точный ответ, сможете ли вы реализовать приближение под ваши ограничения или нет? Это похоже на проблему программирования сейчас. μ+loglogμlogμμ+pμlogμ
Дуглас Заре

Ответы:

10

Распределение Пуассона с большим средним значением приблизительно нормально, но вы должны быть осторожны, если хотите, чтобы хвост был ограничен, а нормальное приближение пропорционально менее точно вблизи хвостов.

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

k=Dexp(μ)μkk!<k=Dexp(μ)μDD!(μD+1)kD=exp(μ)μDD!11μD+1<exp(μ)μD2πD(D/e)D11μD+1=exp(Dμ)(μD)DD+12πD(D+1μ)

Линия 2 строка 3 была связана с формулой Стирлинга. На практике я думаю, что вы хотите решить - p log 2 = log ( bound ) численно с помощью бинарного поиска. Метод Ньютона, начинающийся с начального предположения D = μ + c plog2=log(bound)также должен работать.D=μ+cμ.

Например, при и µ = 1000 численное решение, которое я получаю, составляет 1384,89. Распределение Пуассона со средним значением 1000 принимает значения от 0 до 1384p=100μ=100010000138411/2100.06.0138311/299.59.

Дуглас Заре
источник
1
+1. Другой подход связывает вероятности хвоста Пуассона (справа) с вероятностями хвоста гамма-распределений (слева), которые можно близко (более) оценить в приближении седловой точки.
whuber
От этого далеко до чего-то, ограниченного 64-битной целочисленной арифметикой (без exp, log, sqrt ..), но я над этим поработаю; Спасибо всем!
fgrieu
(+1) Вплоть до вызова приближения Стерлинга (что не имеет значения) это как раз та граница, на которую я (непрозрачно) ссылался в своем комментарии к ОП. (Например, см. Здесь .)
кардинал
2

Yλ

G(x)=2(xlnxλ+λx)  sign(xλ).
Φk0
п(Y<К)Φ(г(К))п(YК),
Φ(г(К-1))п(Y<К)Φ(г(К))
для всех целых К>0 . Кроме того, Φ(G(k+(1/2)))P(Yk)Y < K ) Ф
Φ(G(k1/2))P(Y<k)Φ(G(k))
k>0

Павел Рузанкин
источник
Если бы вы могли выписать ключевое уравнение (при условии, что их всего один или два), это помогло бы, если бы в какой-то момент связь перестала работать.
jbowman