14

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

Improve this question

8 Answers 8

Reset to default
36

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
Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

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.
you had me at "Note that that would generate the interarrival times".. best explanation that I have seen on this
7

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.

Improve this answer

4 Comments

Isn't this is a link to a Poisson distribution, not a Poisson process?
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.
By the way, isn't the mean (scale) of the exponential distribution 1/k, not k?
There are two different conventions for parameterizing the exponential distribution.
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

Improve this answer

1 Comment

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

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.

Improve this answer

2 Comments

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

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

Improve this answer

Comments

1

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

Improve this answer

Comments

0

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://stackoverflow.com/a/15307412/1650437

Improve this answer

Comments

0

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)
Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.