How do I generate a Poisson Process?

Question

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

Aldana

Answer 1

Here's sample code for generating Poisson samples using C++ TR1.

If you want a Poisson process, times between arrivals are exponentially distributed, and exponential values can be generated trivially with the inverse CDF method: -k*log(u) where u is a uniform random variable and k is the mean of the exponential.

Noblewoman

Isn't this is a link to a Poisson distribution, not a Poisson process? – Aldana

Good point. If you want a Poisson process, times between arrivals are exponentially distributed, and exponential values can be generated trivially with the inverse CDF method: -k*log(u) where u is a uniform random variable and k is the mean of the exponential. – Noblewoman

By the way, isn't the mean (scale) of the exponential distribution 1/k, not k? – Amino

There are two different conventions for parameterizing the exponential distribution. – Noblewoman

Answer 2

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

Sinciput

Be sure not to take log of 0 when calculating log(rand()). A common trick is to calculate log(1.0 - rand()) instead, as rand() usually returns a number, which is less than 1. – Detail

you had me at "Note that that would generate the interarrival times".. best explanation that I have seen on this – Elwina

Answer 3

I would be very careful about using the inverse CDF and pumping a uniform random number through it. The problem here is that often the inverse CDF is numerically unstable or the functions to produce it can give undesirable fluctuations near the ends of the interval. For that reason I would recommend something like the rejection method used in "Numerical Recipes in C". See the poidev function given in ch 7.3 of NRC: http://www.nrbook.com/a/bookcpdf/c7-3.pdf

Linnette

The link is to a blank pdf. I have the C++ version (v3 I believe) and it doesn't even include a Poisson deviate. But, my understanding is that the deviates are samples from the distribution which would leave me where I started. – Aldana

Answer 4

Generating arrival times via Poisson Process does not mean using a Poisson distribution. It is done by creating an exponential distribution based on the Poisson arrival rate lamda.

In short, you need to generate an exponential distribution with an average = 1/lamda, see the following example:

#include <iostream>
#include <iterator>
#include <random>

int
main ()
{
 // seed the RNG
 std::random_device rd; // uniformly-distributed integer random number generator
 std::mt19937 rng (rd ()); // mt19937: Pseudo-random number generation

 double averageArrival = 15;
 double lamda = 1 / averageArrival;
 std::exponential_distribution<double> exp (lamda);

double sumArrivalTimes=0;
double newArrivalTime;


 for (int i = 0; i < 10; ++i)
  {
   newArrivalTime=  exp.operator() (rng); // generates the next random number in the distribution 
   sumArrivalTimes  = sumArrivalTimes + newArrivalTime;  
   std::cout << "newArrivalTime:  " << newArrivalTime  << "    ,sumArrivalTimes:  " << sumArrivalTimes << std::endl;  
  }

}

The result of running this code:

newArrivalTime:  21.6419    ,sumArrivalTimes:  21.6419
newArrivalTime:  1.64205    ,sumArrivalTimes:  23.2839
newArrivalTime:  8.35292    ,sumArrivalTimes:  31.6368
newArrivalTime:  1.82962    ,sumArrivalTimes:  33.4665
newArrivalTime:  34.7628    ,sumArrivalTimes:  68.2292
newArrivalTime:  26.0752    ,sumArrivalTimes:  94.3045
newArrivalTime:  63.4728    ,sumArrivalTimes:  157.777
newArrivalTime:  3.22149    ,sumArrivalTimes:  160.999
newArrivalTime:  1.64637    ,sumArrivalTimes:  162.645
newArrivalTime:  13.8235    ,sumArrivalTimes:  176.469

so, based on your experiment you can either use: newArrivalTime or sumArrivalTimes.

ref: http://www.math.wsu.edu/faculty/genz/416/lect/l05-45.pdf

Tko

Answer 5

If you are using python, you can use random.expovariate(rate) to generate arrival times at rate events per time interval

Venation

Answer 6

In order to pick a sample from a distribution, you need to compute the inverse cumulative distribution function (CDF). You first pick a random number uniformly on the real interval [0, 1], and then take the inverse CDF of that value.

Drugi

Thanks - my thoughts were just that, but I felt strange since the cdf generates N(t) but not t. And, I didn't know the correct way to sample from N(t) to get a Poisson process. What happens if I don't sample uniformly from the CDF (will it still be a Poisson process)? – Aldana

Inverting the CDF of a Poisson is not easy or efficient. For a more efficient approach, see Corwin's link or see my answer about how to use C++ TR1. – Noblewoman

Answer 7

In python, you can try below code.

If you want to generate 20 random readings in 60 seconds. ie (20 is the lambda)

 def poisson_job_generator():
    rateParameter = 1.0/float(60/20) 
    while True:
        sl = random.expovariate(rateParameter)

Arlin

Answer 8

The discussion here has all the details about using inverse sampling to generate inter-arrivals, which is usually what people want to do for games.

https://mcmap.net/q/827096/-simulating-poisson-waiting-times

Staggard