如何比这更快地对 Python 中的非齐次泊松过程进行采样？
How to sample inhomogeneous Poisson processes in Python faster than this?
我正在对速率不固定的毫秒时间尺度的泊松过程进行采样。我通过检查每个大小为 delta 的间隔是否存在基于该间隔中的平均速率的事件来离散化采样过程。由于我使用的是 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)))
这是一个运行速度提高大约 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
旁注：不过，这并不是模拟非齐次泊松过程的最佳方法。来自 Wikipedia:
An inhomogeneous Poisson process with intensity function λ(t) can be simulated by rejection sampling from a homogeneous Poisson process with fixed rate λ: choose a sufficiently large λ so that λ(t) = λ p(t) and simulate a Poisson process with rate parameter λ. Accept an event from the Poisson simulation at time t with probability p(t).
我正在对速率不固定的毫秒时间尺度的泊松过程进行采样。我通过检查每个大小为 delta 的间隔是否存在基于该间隔中的平均速率的事件来离散化采样过程。由于我使用的是 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)))
这是一个运行速度提高大约 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
旁注：不过，这并不是模拟非齐次泊松过程的最佳方法。来自 Wikipedia:
An inhomogeneous Poisson process with intensity function λ(t) can be simulated by rejection sampling from a homogeneous Poisson process with fixed rate λ: choose a sufficiently large λ so that λ(t) = λ p(t) and simulate a Poisson process with rate parameter λ. Accept an event from the Poisson simulation at time t with probability p(t).