如何以比这更快的速度在Python中采样非均匀泊松过程？

Question

我正在以不固定速率的毫秒级采样泊松过程。我通过基于该间隔的平均速率检查每个大小增量的间隔中是否存在事件来离散化采样过程。由于我正在使用Python，因此它的运行速度比我希望的要慢一些。我当前使用的代码如下：

import numpy
def generate_times(rate_function,max_t,delta):
    times = []
    for t in numpy.arange(delta,max_t,delta):
        avg_rate = (rate_function(t)+rate_function(t+delta))/2.0
        if numpy.random.binomial(1,1-math.exp(-avg_rate*delta/1000.0))>0:
            times.extend([t])
    return times

速率函数可以是任意的，考虑到速率函数，我不是在寻找封闭形式的解决方案。

如果您想使用一些参数，可以尝试：

max_t = 1000.0
delta = 0.1
high_rate = 100.0
low_rate = 0.0
phase_length = 25.0
rate_function = (lambda x: low_rate + (high_rate-low_rate)*0.5*(1+math.sin(2*math.pi*x/phase_length)))

Answer 1

这是一个运行速度快约75倍并实现相同功能的版本：

def generate_times_opt(rate_function,max_t,delta):
    t = numpy.arange(delta,max_t, delta)
    avg_rate = (rate_function(t) + rate_function(t + delta)) / 2.0
    avg_prob = 1 - numpy.exp(-avg_rate * delta / 1000.0)
    rand_throws = numpy.random.uniform(size=t.shape[0])

    return t[avg_prob >= rand_throws].tolist()

输出：

11:59:07 [70]: %timeit generate_times(rate_function, max_t, delta)
10 loops, best of 3: 75 ms per loop

11:59:23 [71]: %timeit generate_times_opt(rate_function, max_t, delta)
1000 loops, best of 3: 1.15 ms per loop

旁注：不过，这不是模拟非均匀泊松过程的最佳方法。从维基百科：

强度函数为ω（t）的非均匀泊松过程可以通过以固定速率ω从均匀泊松过程中进行拒绝采样来模拟：选择一个足够大的σ。这样？（t）=？p（t）并用速率参数？模拟泊松过程。在时间t以概率p（t）接受来自Poisson模拟的事件。