Как моделировать из распределения Пуассона, используя моделирование из экспоненциального распределения

Меня просят реализовать алгоритм моделирования распределения Пуассона (лямбда) с использованием моделирования экспоненциального распределения.

Мне была задана следующая плотность: P(X = k) = P(X1 + · · · + Xk ≤ 1 ‹ X1 + · · · + Xk+1), для k = 1, 2, . . . . P(X = k) — пуассон с лямбдой, а Xi — экспоненциальное распределение.

Я написал код для имитации экспоненциального распределения, но понятия не имею, как имитировать пуассон. Может ли кто-нибудь помочь мне в этом? Спасибо миллион.

Мой код:

n<-c(1:k)
  u<-runif(k)
  x<--log(1-u)/lambda

r poisson simulation
person Amy_777    schedule 23.10.2017    source источник
comment
используйте функцию rpois для генерации случайных чисел для распределения Пуассона с параметром лямбда.   -  person Samuel    schedule 23.10.2017

Ответы (2)


arrow_upward
3
arrow_downward

Я работаю исходя из предположения, что вы (или ваш инструктор) хотите сделать это из первых принципов, а не просто вызывать встроенный генератор Пуассона. Алгоритм довольно прост. Вы подсчитываете, сколько экспонент вы можете сгенерировать с указанной скоростью, пока их сумма не превысит 1.

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

count <- 0
sum <- 0
repeat {
  generate x ~ exp(lambda)
  sum <- sum + x
  if sum > 1
    break
  else
    count <- count + 1
}

Значение count после того, как вы break вышли из цикла, является результатом Пуассона для этого испытания. Если вы оберните это как функцию, верните count, а не breaking из цикла.

Вы можете улучшить это вычислительно несколькими способами. Во-первых, обратите внимание, что член 1-U для генерации экспонент имеет равномерное распределение и может быть заменен только U. Более значительное улучшение достигается при записи оценки как максимизировать i s.t. SUM(-log(Ui) / rate) <= 1, значит SUM(log(Ui)) >= -rate.

Теперь возведите в степень обе стороны и упростите, чтобы получить

PRODUCT(Ui) >= Exp(-rate).

Правая часть этого числа постоянна и может быть вычислена заранее, уменьшая объем работы от k+1 вычислений журнала и дополнений до одного возведения в степень и k+1 умножений:

count <- 0
product <- 1
threshold = Exp(-lambda)
repeat {
  generate u ~ Uniform(0,1)
  product <- product * u
  if product < threshold
    break
  else
    count <- count + 1
}

Предполагая, что вы выполняете замену U на 1-U для обеих реализаций, они алгебраически равны и дадут идентичные ответы с точностью до арифметики с плавающей запятой для заданного набора U.

person pjs    schedule 23.10.2017

arrow_upward
2
arrow_downward

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

Шаг 1. Создайте (большую) выборку из экспоненциального распределения и создайте вектор кумулятивных сумм. k-й элемент этого вектора есть время ожидания k-го пуассоновского прихода

Шаг 2. Измерьте, сколько прибытий мы наблюдаем за единичный интервал времени.

Шаг 3. Многократно повторите шаги 1 и 2 и соберите результаты в вектор.

Это будет ваша выборка из распределения Пуассона с правильным параметром скорости.

Код:

lambda=20 # for example
out=sapply(1:100000, function(i){
   u<-runif(100)
   x<--log(1-u)/lambda
   y=cumsum(x)
   length(which(y<=1))
})

Затем вы можете проверить достоверность встроенной функции с помощью теста Колмогорова-Смирнова:

ks.test(out, rpois(100000, lambda))
person ags29    schedule 23.10.2017