How to sample inhomogeneous Poisson processes in Python faster than this?

Question

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

Responsiveness

Are you asking how to optimize generate_times? Can it return an ndarray or must it be a list? – Engle

It does not need to be fully optimised, any factor greater than 5 in the runtime would make a huge difference. And the output must be a list (you can convert an ndarray into a list in the end if that's faster... :P). – Responsiveness

generate_times seems incomplete, it never adds anything to tunes – Engle

Thanks, fixed it now. – Responsiveness

Have you looked at other SO answers regarding vectorizing calcs on ndarrays? – Engle

If the rate is not constant why do you call it a Poisson process? I thought that a constant arrival rate was part of the definition of a Poisson process. – Concubine

@JohnColeman en.wikipedia.org/wiki/Inhomogeneous_Poisson_process – Impute

@Impute I see. It is a generalization of the classic Poisson process. Thanks for the link. – Concubine

Answer 1

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

Impute

Awesome, thanks a lot, the alternative from the wikipedia page also makes a lot of sense. It is weak in respect to having one maximum for the rate function and other values close to 0. – Responsiveness