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

Are you asking how to _optimize_ ```generate_times```? Can it return an ndarray or must it be a list? — wwii, Sep 22 '15 at 09:34
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). — Haffi112, Sep 22 '15 at 09:38
```generate_times``` seems incomplete, it never adds anything to ```tunes``` — wwii, Sep 22 '15 at 09:38
Have you looked at other SO answers regarding _vectorizing_ calcs on ndarrays? — wwii, Sep 22 '15 at 09:58
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. — John Coleman, Sep 22 '15 at 10:12
@JohnColeman https://en.wikipedia.org/wiki/Inhomogeneous_Poisson_process — musically_ut, Sep 22 '15 at 10:50
@musically_ut I see. It is a generalization of the classic Poisson process. Thanks for the link. — John Coleman, Sep 22 '15 at 11:03

musically_ut · Accepted Answer · 2015-09-22T10:52:57.433

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

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. — Haffi112, Sep 22 '15 at 11:01

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

1 Answers1