3

我有 N 个进程，每个进程都有不同的泊松率。我想模拟所有 N 个进程的到达时间。如果 N =1 我可以这样做

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

但是，如果我有N = 5 费率清单

rates =  [10,1,15,4,2]

我希望循环以正确的顺序输出所有 N 个进程的到达时间。也就是说，我仍然希望每行只有两个数字（进程的 ID 和到达时间），但按到达时间全局排序。

我可以制作 N 个列表并在之后合并它们，但我希望首先以正确的顺序输出到达时间。

更新。一个问题是，如果您只是从每个流程中采样固定数量的到达，您只能从高速率流程中获得早期时间。所以我认为我需要从每个过程的固定时间间隔进行采样，因此样本数量会根据速率而变化。

4

2 回答 2

3

如果我对您的理解正确：

import random
import itertools

def arrivalGenerator(rate):
    t = 0
    while True:
        t += random.expovariate(rate)
        yield t

rates = [10, 1, 15, 4, 2]
t = [(i, 0) for i in range(0, len(rates))]
arrivals = []
for i in range(len(rates)):
    t = 0
    generator = arrivalGenerator(rates[i])
    arrivals += [(i, arrival) \
                 for arrival in itertools.takewhile(lambda t: t < 100, generator)]


sorted_arrivals = sorted(arrivals, key=lambda x: x[1])
for arrival in sorted_arrivals:
    print arrival[0], arrival[1]

请注意，您的初始逻辑是为每个进程生成固定数量的到达。您真正想要的是一个特定的时间窗口，并继续为给定进程生成，直到您超过该时间窗口。

于 2013-07-05T14:52:54.043 回答
1

在http://www.columbia.edu/~ks20/4703-Sigman/4703-07-Notes-PP-NSPP.pdf之后，我认为有一个更有效的答案。

你大致做：

total_rate = sum(rates)
probabilities = [ r/total_rate for r in rates ]

arrivals = []
t = 0
while t < T:
   t += random.expovariate(total_rate)
   i = weighted_random(probabilities)
   arrivals += (i, t)

这种方法消除了为大量不同的到达过程保持协程状态的需要。只有一个“净”到达过程。分配将是相同的。

请注意，我没有在上面给出 weighted_random 的实现，但我认为我的意图很明确。它留给读者作为练习;-) - 或参见例如http://eli.thegreenplace.net/2010/01/22/weighted-random-generation-in-python。

你也可以这样做：

arrivals = []
t = 0
while t < T:
    dt_list = [ random.expovariate(r) for r in rates ]
    (dt,i) = min((tau,i) for i,tau in enumerate(dt_list))
    t += dt
    arrivals += (i, t)

即，您实际上确实为所有进程生成了单独的到达间隔时间，但您不需要“记住”进程的状态。请注意，具有速率 r1 和 r2 的两个独立指数分布随机变量的最小值本身以速率 r1+r2 呈指数分布（根据http://ocw.mit.edu/courses/mathematics/18-440-probability-and- random-variables-spring-2011/lecture-notes/MIT18_440S11_Lecture20.pdf），所以这实际上与前面的代码片段非常相似。

在我在这里给出的两种方法中，我认为第一种更好：

  • 第一个是 O( len(arrivals) * log(len(rates)) ) 而第二个是 O( len(arrivals) * len(rates) )
  • 第一个每次到达需要来自底层生成器的 2 个随机数，而第二个每次到达需要 len(rates) 个随机数。
  • 第一个需要每次到达对 log 进行 1 次评估（我假设这是生成指数随机变量的方式），而第二个需要每次到达对 log 进行 O(len(rates)) 评估。

此外，对上述所有 Python 语法持保留态度（我没有运行它，而且我对 Python 很生疏），如果你愿意，可以消除临时列表。这实际上是“伪代码”；对于快速的蒙特卡罗模拟，无论如何您可能都会使用 C++（和/或 CUDA）。

我知道你可能已经过了需要这个答案的地步，但我希望它对找到这篇文章的其他人有所帮助。

于 2015-01-15T16:19:05.513 回答