Recently I got it into my head that I wanted to be able to generate random numbers from a Poisson distribution. I looked around and found a number of topics on how to generate a table of probabilities from a given distribution and then use it to look up correspondence from uniform random numbers. This probably works wery well but I wanted to have something faster that didn't require such an amount of precomputing for every distribution. The problem, if I've got it correctly is that it's not all that easy to invert the Poisson cdf in order to map values - and apart from that the distribution is descrete so the mapping wouldn't even be one-to-one. Then I found that some people recommended using exponential distributions instead of Poisson if possible, since it is trivial to map uniformly distributed random values to exponentially distributed ones. Now this is nice but for the fact that when comparing Poisson to exponential distributions (ignoring the fact that one is discrete while the other is continous) the exponential has, regardless of expectation, a non-negligible probability that the value is close to zero, whereas the Poisson is more similar to a normal distribution with larger expectation. Put in other words the variance of the exponential distribution is the square of the expectation while with Poisson they are equal. This is one of the properties of Poisson that I find useful and I'm not ready to trade it off for simpler computation. Finally I stumbled over a rumor (without reference) that using normally distributed values (generated with Box-Muller or something) and transforming these to exponentially distributed ones I would get something similar to values of a Poisson distribution. Is there any truth to this? I must admit that intuitively it seems possible as long as one keeps a close eye at expectation and variance along the way. But I can't back it up with any theory. Would this be a good approximation, and in such a case, within what boundaries? What methods do you use to generate non-uniform random numbers? Thanks for reading this far.

The continuous exponential probability distribution is related to the discrete Poisson distribution. The Poisson distribution provides an appropriate description of the number of occurrences per interval. The exponential distribution provides a description of the length of the interval between occurrences.

There is a function called poidev() in Numerical Recipes which generates random (integer) values that form a poisson distribution. You can find it here (Chapter 7.3). But in case to use it you need this (6.1 for gammln()) and this (7.1 for ran1()).

I'm currently using it in a project, and despite a little miss in the distribution it gives good results.

Hope that helps.

HellRaiZer

Quote:Original post by y2kiah
The continuous exponential probability distribution is related to the discrete Poisson distribution. The Poisson distribution provides an appropriate description of the number of occurrences per interval. The exponential distribution provides a description of the length of the interval between occurrences.

I spent some time last evening contemplating this. I wanted to verify that it held in practice so I wrote a little function that
1. generates a uniform random number via a library routine,
2. transforms the uniform value to an exponentially distributed one, and
3. adds the exponentially distributed value to a sum if it doesent exceed a given interval. (Unit interval in this case)
Then it repeats the process until the interval is exceeded and returns the number of iterations.

I found by generating lots of these values and presenting them in a histogram, that the values generated fits very well with the poisson distribution. The function that I used can be seen below. Now my followup question is, should this method be used in practice? I understand that it is inefficient for large mu (even in risk of stalling due to underflow in the addition to sum). But if the method is extended to use a normal distribution approximation for mu larger than say 20, wouldn't that mend the issue?

Anyway, thanks to y2kiah for the conceptual input.

Trap and HellRaiZer: I've glanced over the linked pages both of you posted but I'll have to get back to you once I've read them more carefully.

The code I used: (in MATLAB syntax if you wonder)