views:

1182

answers:

3

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 lambda 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 lambda) I could use Poisson(t2 x lambda)-Poisson(t1 x lambda) 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).

+1  A: 

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.

Adam Rosenfield
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)?
bias
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.
John D. Cook
+3  A: 

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

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.
bias
+3  A: 

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

John D. Cook
Isn't this is a link to a Poisson distribution, not a Poisson process?
bias
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.
John D. Cook