13

原始问题：

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

澄清：

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

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

我可以通过使间隔[t2,t1]足够小来做到这一点，以便每个间隔只有一次到达（发生为|t2-t1| -> 0）。

4

8 回答 8

31

如果您有一个速率参数为 L 的泊松过程（这意味着，从长远来看，每秒有 L 个到达），那么到达间隔时间呈指数分布，均值为 1/L。所以 PDF 为 f(t) = -L*exp(-Lt)，CDF 为 F(t) = Prob(T < t) = 1 - exp(-Lt)。所以你的问题变成了：如何生成一个分布为 F(t) = 1 - \exp(-Lt) 的随机数 t？

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

-log(rand()) / L

由于 python 提供了一个函数来生成指数分布的随机数，您可以模拟泊松过程中的前 10 个事件，平均每秒到达 15 次，如下所示：

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

请注意，这将生成 *inter*arrival 时间。如果您想要到达时间，则必须像这样继续向前移动时间变量：

import random
t= 0
for i in range(1,10):
   t+= random.expovariate(15)
   print t
于 2012-04-20T17:35:44.200 回答
6

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

如果您想要泊松过程，到达之间的时间呈指数分布，并且可以使用逆 CDF 方法轻松生成指数值：-k*log(u) 其中 u 是均匀随机变量，k 是指数的平均值。

于 2009-07-21T00:40:19.057 回答
3

我会非常小心地使用逆 CDF 并通过它抽取一个统一的随机数。这里的问题是，逆 CDF 通常在数值上是不稳定的，或者产生它的函数会在区间末端附近产生不希望的波动。出于这个原因，我会推荐类似“C 中的数字食谱”中使用的拒绝方法。请参阅 NRC 第 7.3 章中给出的 poidev 函数：http ://www.nrbook.com/a/bookcpdf/c7-3.pdf

于 2009-07-20T21:15:38.580 回答
1

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

于 2009-07-20T19:56:12.420 回答
1

如果您使用的是 python，则可以使用 random.expovariate(rate) 以每个时间间隔的速率事件生成到达时间

于 2011-04-15T04:55:07.547 回答
0

这里的讨论包含有关使用反向采样生成到达间隔的所有细节，这通常是人们想要为游戏做的事情。

https://stackoverflow.com/a/15307412/1650437

于 2013-03-09T05:15:20.190 回答
0

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

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

 def poisson_job_generator():
    rateParameter = 1.0/float(60/20) 
    while True:
        sl = random.expovariate(rateParameter)
于 2017-05-11T21:49:51.910 回答
0

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

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

#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

于 2017-12-07T16:16:58.377 回答