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

I'm sampling a Poisson process at a millisecond time scale where the rate is not fixed. I discretise the sampling process by checking in each interval of size delta whether there is an event there or not based on the average rate in that interval. Since I'm using Python it's running a bit slower than I would hope it to be. The code I'm currently using is the following:

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

The rate function can be arbitrary, I'm not looking for a closed form solution given a rate function.

If you want some parameters to play with you can try:

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)))

解决方案

Here's a version which runs about 75x faster and implements the same function:

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()

Output:

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

Sidenote: This is not the best way to simulate an Inhomogenous Poisson Process, though. From 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).