如何绘制带指数核的泊松过程

Question

我想在一个特定的时间范围内模拟一个泊松过程的时间，这个过程有一个指数核。我有一段代码可以实现这个功能，但写得非常糟糕。

# Attempt to simulate a Poisson process.
import random
import matplotlib.pyplot as plt
t = 0
rate = 5
timeinterval = 100
t = int(100*random.expovariate(rate))
times=[]
while (t < timeinterval):
    times.append(t)
    print t
    t+= int(100*random.expovariate(rate))

print times
s = [0]*100
a=0.5
for i in xrange(len(s)):
    s[i] = int(i-1 in times) + (1-a) * s[i-1]
plt.plot(s)
plt.show()

我该如何避免这段代码中的一些不太优雅的写法呢？

举个例子，我为了能画出数据图，乘了100并使用了整数类型，但其实我不想这样做。

另外，我希望时间间隔是1而不是100（已经把代码中所有乘以100的部分去掉了），但我不知道该怎么画出结果。

Answer 1

这些小问题的主要原因是，你通过采样尖峰出现的时间来解决问题，但接着你又想把它们强行放到一个网格上，因为你希望它们逐渐减弱。

所以我建议从一个网格开始，计算尖峰发生的概率。这样你就可以同时让信号减弱：

import random
import matplotlib.pyplot as plt
from numpy import linspace, exp

random.seed(123)
t0 = 0
t1 = 1
times = linspace(t0, t1, 1000)
dt = times[1] - times[0]
result = []
s = 0
decay = 0.01  # decay to 1/e after this time                                                       
lambd = 0.2   # expect roughly one spike each timeinterval lambd                                   
n = 0
for t in times:
    if random.random() > exp(-dt/lambd): # actually random.random() < 1 - exp(-dt/lambd)           
        s += 1.0
        n += 1
    s = s * exp(-dt/decay)
    result.append(s)

print "number of observed spikes=%i expected mean=%f" % (n, (t1-t0)/lambd)
plt.figure(figsize=(6,4))
plt.plot(times, result)
plt.savefig("process.png")
plt.show()

这是种子 123 的一个示例输出：

请注意，无论你改变 linspace 中最后一个参数的值，得到的结果都是类似的（虽然由于随机过程，结果不完全相同），这正是我们所期望的。

补充：你也可以使用 expovariate 来实现，这更接近你最初的方法。

rate = 1.0/lambd
tspike = t0 + random.expovariate(rate)
for t in times:
    s = s * exp(-dt/decay)
    if t > tspike:
        s += 1.0
        n += 1
        tspike = t + random.expovariate(rate)
    result.append(s)

在这里，我还交换了减弱和尖峰生成的顺序，所以它们从 1 开始。