How do I generate a Poisson Process?

bias picture bias · Jul 20, 2009 · Viewed 28.1k times · Source

Original Question:

I want to generate a Poisson process. If the number of arrivals by time t is N(t) and I have a Poisson distribution with parameter λ how do I generate N(t)? How would I do this in C++?

Clarification:

I originally wanted to generate the process using a Poisson distribution. But, I was confused about what parameter from the process I needed; I thought I could use N(t) but that tells me how many arrivals have occurred on the interval (0,t] which wasn't what I wanted. So, then I thought I could use N(t2)-N(t1) to get the number of arrivals on the interval [t1,t2]. Since N(t)~Poisson(t x λ) I could use Poisson(t2 x λ)-Poisson(t1 x λ) but I don't want the number of arrivals in an interval.

Rather, I want to generate the explicit times that arrivals occur at.

I could do this by making the interval [t2,t1] sufficiently small so that each interval has only one arrival (which occurs as |t2-t1| -> 0).

Answer

Chris Marshall picture Chris Marshall · Apr 20, 2012

If you have a Poisson process with rate parameter L (meaning that, long term, there are L arrivals per second), then the inter-arrival times are exponentially distributed with mean 1/L. So the PDF is f(t) = -L*exp(-Lt), and the CDF is F(t) = Prob(T < t) = 1 - exp(-Lt). So your problem changes to: how to I generate a random number t with distribution F(t) = 1 - \exp(-Lt)?

Assuming the language you are using has a function (let's call it rand()) to generate random numbers uniformly distributed between 0 and 1, the inverse CDF technique reduces to calculating:

-log(rand()) / L

As python provides a function to generate exponentially distributed random numbers, you could simulate the first 10 events in a poisson process with an averate rate of 15 arrivals per second like this:

import random
for i in range(1,10):
   print random.expovariate(15)

Note that that would generate the *inter*arrival times. If you wanted the arrival times, you would have to keep moving a time variable forward like this:

import random
t= 0
for i in range(1,10):
   t+= random.expovariate(15)
   print t