我该如何生成泊松过程？

Question

原始问题：

我想生成泊松过程。如果时间t内到达的数量为N(t)，并且我有一个参数为λ的泊松分布，那么如何生成N(t)？我该如何在C++中实现这一点？

澄清：

我最初想使用泊松分布来生成过程。但是，我对需要哪个过程参数感到困惑；我以为我可以使用N(t)，但这告诉我已经发生了多少到达在区间(0,t]，这不是我想要的。所以，然后我想我可以使用N(t2)-N(t1)来获得区间[t1,t2]上到达的数量。由于N(t)~Poisson(t x λ)，我可以使用Poisson(t2 x λ)-Poisson(t1 x λ)，但我不想知道区间内到达的数量。

相反，我想生成到达发生的明确时间。

我可以通过将区间[t2,t1]足够小，以便每个区间只有一个到达（出现在|t2-t1| -> 0）来实现这一点。

- bias

Answer 1

如果您有一个速率参数为L的泊松过程（长期来看，每秒钟有L个到达），那么到达间隔时间呈指数分布，平均值为1/L。因此，概率密度函数为f(t)=-L*exp(-Lt)，累积分布函数为F(t)=Prob(T

假设您使用的语言有一个函数（我们称之为rand()）可以生成在0和1之间均匀分布的随机数，则反函数CDF技术可简化为计算：

-log(rand()) / L

由于Python提供了生成指数分布随机数的函数，您可以像这样模拟平均到达率为每秒15个的泊松过程中的前10个事件：

import random
for i in range(1,10):
   print random.expovariate(15)

请注意，这将生成*间隔时间*。如果您想要到达时间，您需要像这样不断地向前移动时间变量：

import random
t= 0
for i in range(1,10):
   t+= random.expovariate(15)
   print t

- Chris Marshall

---

在计算log(rand())时，一定不要取0的对数。一个常见的技巧是计算log(1.0 - rand())，因为rand()通常返回小于1的数字。 - fdermishin

你说的“注意这将生成间隔时间”一语击中了我，这是我见过的最好的解释。 - ihightower

Answer 2

这是使用 C++ TR1 生成泊松样本的示例代码。

如果你想要一个泊松过程，到达时间之间的时间服从指数分布，可以通过反函数法轻松生成指数值：-k*log(u)，其中u为均匀随机变量，k为指数的平均值。

- John D. Cook

---

这不是一个泊松过程的链接吗？为什么是泊松分布的链接？ - bias

好的观点。如果你想要一个泊松过程，到达时间间隔是指数分布的，而指数值可以通过反函数法轻松生成：-k*log(u)，其中u是均匀随机变量，k是指数分布的平均值。 - John D. Cook

顺便问一下，指数分布的均值（尺度）不是1/k而是k吗？ - Peter O.

指数分布有两种不同的参数化约定。 - John D. Cook

Answer 3

关于使用逆CDF并将均匀随机数通过它的问题，我建议非常谨慎。这里的问题是逆CDF经常是数值不稳定的，或者产生它的函数可能在区间的两端产生不良波动。因此，我建议使用“Numerical Recipes in C”中使用的拒绝方法。请参见NRC第7.3章中给出的poidev函数：http://www.nrbook.com/a/bookcpdf/c7-3.pdf

- Corwin Joy

---

链接是一个空白的PDF文件。我有C++版本（我相信是v3），甚至没有包括泊松偏差。但是，我的理解是这些偏差是从分布中抽样得到的，这会让我回到起点。 - bias

Answer 4

如果您正在使用Python，可以使用random.expovariate(rate)函数来生成每个时间间隔内以rate事件为速率的到达时间。

- Steve Brown

Answer 5

为了从分布中选取一个样本，您需要计算反累积分布函数（CDF）。首先，在实数区间[0,1]上均匀地选择一个随机数，然后取该值的反CDF。

- Adam Rosenfield

---

谢谢 - 我的想法就是这样，但因为累积分布函数生成的是 N(t) 而不是 t，所以感觉有点奇怪。同时，我也不知道从 N(t) 中采样得到泊松过程的正确方法。如果我没有均匀地从累积分布函数中采样，会发生什么（它仍然是泊松过程吗）？ - bias

将泊松分布的累积分布函数反转并不容易或高效。要获得更高效的方法，请参见Corwin的链接，或查看我的答案以了解如何使用C++ TR1。 - John D. Cook

Answer 6

通过泊松过程生成到达时间并不意味着使用泊松分布。这是通过基于泊松到达率λ创建指数分布来完成的。

简而言之，您需要生成一个平均值为1/λ的指数分布，请参见以下示例：

#include <iostream>
#include <iterator>
#include <random>

int
main ()
{
 // seed the RNG
 std::random_device rd; // uniformly-distributed integer random number generator
 std::mt19937 rng (rd ()); // mt19937: Pseudo-random number generation

 double averageArrival = 15;
 double lamda = 1 / averageArrival;
 std::exponential_distribution<double> exp (lamda);

double sumArrivalTimes=0;
double newArrivalTime;


 for (int i = 0; i < 10; ++i)
  {
   newArrivalTime=  exp.operator() (rng); // generates the next random number in the distribution 
   sumArrivalTimes  = sumArrivalTimes + newArrivalTime;  
   std::cout << "newArrivalTime:  " << newArrivalTime  << "    ,sumArrivalTimes:  " << sumArrivalTimes << std::endl;  
  }

}

运行此代码的结果：

newArrivalTime:  21.6419    ,sumArrivalTimes:  21.6419
newArrivalTime:  1.64205    ,sumArrivalTimes:  23.2839
newArrivalTime:  8.35292    ,sumArrivalTimes:  31.6368
newArrivalTime:  1.82962    ,sumArrivalTimes:  33.4665
newArrivalTime:  34.7628    ,sumArrivalTimes:  68.2292
newArrivalTime:  26.0752    ,sumArrivalTimes:  94.3045
newArrivalTime:  63.4728    ,sumArrivalTimes:  157.777
newArrivalTime:  3.22149    ,sumArrivalTimes:  160.999
newArrivalTime:  1.64637    ,sumArrivalTimes:  162.645
newArrivalTime:  13.8235    ,sumArrivalTimes:  176.469

所以，根据你的实验，你可以使用 newArrivalTime 或 sumArrivalTimes。

参考：http://www.math.wsu.edu/faculty/genz/416/lect/l05-45.pdf

- Moe

Answer 7

这里的讨论详细介绍了使用反向抽样生成间隔时间的所有细节，这通常是游戏开发者想要做的。

https://dev59.com/0ljUa4cB1Zd3GeqPUMLV#15307412

- jdbertron

Answer 8

在Python中，您可以尝试以下代码。

如果您想在60秒内生成20个随机读数。即（20是lambda）

 def poisson_job_generator():
    rateParameter = 1.0/float(60/20) 
    while True:
        sl = random.expovariate(rateParameter)

- josepainumkal