如何用指数核绘制泊松过程答案

【问题标题】：How to plot a Poisson process with an exponential kernel如何用指数核绘制泊松过程
【发布时间】：2014-08-29 16:30:51
【问题描述】：

我想在给定时间窗口内模拟具有指数内核的泊松过程的时间。我有以下代码，它有效但很糟糕。

# 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 并使用 int 来绘制数据。理想情况下我不会那样做。

在相关说明中，我希望时间间隔为 1 而不是 100（已经摆脱了代码的所有 *100 部分），但我看不到如何绘制结果。

【问题讨论】：

你为什么要绘制衰变图？您可以在特定时间 t 使用 vlines。那么也没有理由将时间转换为整数网格。
@user333700 我想可视化延迟。最终，我还将在图表中添加其他内容，看看它们如何与延迟相交。
如果您需要延迟，那么在时间点网格上工作是最简单的。但是，您可以使用任何规则间隔的网格来代替整数网格。例如，如果您将times 舍入到小数点后两位，那么您可以在单位间隔 0、0.01、...、1 内构建一个包含 100 个时间点的网格。

标签： python math matplotlib scipy statsmodels

【解决方案1】：

kludges 的主要原因是您通过对尖峰的时间进行采样来解决问题，但是您希望将它们强制放在网格上，因为您希望它们衰减。

因此，我建议从网格开始并计算出现尖峰的概率。然后可以同时衰减信号：

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 中的最后一个参数，无论分辨率如何dt，您都会得到相似（由于随机过程而不同）的结果，就像它应该的那样。

编辑：您也可以使用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开始。

【讨论】：