如何使用指数分布模拟来模拟泊松分布

Question

我被要求实现一种算法，使用指数分布的模拟来模拟泊松（lambda）分布。

\n\n

我得到了以下密度：\nP(X = k) = P(X1 + \xc2\xb7 \xc2\xb7 \xc2\xb7 + Xk \xe2\x89\xa4 1 < X1 + \xc2\xb7 \xc2\ xb7 \xc2\xb7 + Xk+1)，对于 k = 1, 2, . 。。.\nP(X = k) 是具有 lambda 的泊松分布，Xi 是指数分布。

\n\n

我编写了代码来模拟指数分布，但不知道如何模拟泊松分布。有人可以帮我解决这个问题吗？谢谢万。

\n\n

我的代码：

\n\n

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

\n

Answer 1

我正在假设您（或您的讲师）想要从第一原理开始执行此操作，而不仅仅是调用内置的泊松生成器。该算法非常简单。您计算以指定的速率可以生成多少个指数，直到它们的总和超过 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而不是break从循环中进行 ing。

您可以通过多种方式在计算上改进这一点。首先要注意1-U生成指数的项具有均匀分布，并且可以仅替换为U。i通过将评估写为最大化st可以获得更显着的改进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
}

假设您对两种实现进行Ufor1-U替换，它们在代数上是相等的，并且对于给定的一组 ' ，将在浮点算术精度内产生相同的答案U。