问如何利用指数分布模拟泊松分布

EN

我的工作假设你(或者你的导师)想要从基本原则上来做这件事，而不仅仅是调用内置的Poisson生成器。算法非常简单。您可以用指定的速率计算生成多少指数，直到它们的总和超过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的浮点算法的精度范围内得到相同的答案。

按照上述建议，您可以使用rpois生成泊松变量。但是，我对这个问题的理解是，您希望从基本原则出发，而不是使用内置函数。要做到这一点，您需要使用泊松到达的属性，说明到达时间是指数分布的。因此，我们的工作如下：

步骤1:从指数分布生成一个(大)样本，并创建累积和向量。这个向量的第k个入口是第k个泊松到达的等待时间。

步骤2:测量我们在单位时间间隔内看到的到达量

Step3:重复步骤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))
})

然后，您可以通过Kolmogorov测试测试内建函数的有效性：

ks.test(out, rpois(100000, lambda))