How is the poisson distribution a distribution? It seems more like a formula

Question

I just watched this video: https://www.youtube.com/watch?v=Fk02TW6reiA

It shows a formula to calculate an answer for the following problem:

• There are 2 customers expected every 3 minutes in a store
• Therefore there are 6 customers expected every 9 minutes
• What is the likelihood of there being 4 or less in the store in 9 minutes?
• Answer is: P(0;6)+...+P(4;6) which is about 0.28

This makes sense and is well described. However, numpy treats the poisson distribution basically like a random number generator: http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.poisson.html

We can specify lambda as say 5 and how many numbers are desired (the second argument) and get a big list of integers:

>>> import numpy as np
>>> s = np.random.poisson(5, 10000)
>>> s
array([2, 4, 4, ..., 3, 4, 3])
>>> len(s)
10000

These seem like two totally different things. How do you get from using the Poisson formula to calculate the possibility of a certain number of events in a timeframe, to a list of seemingly random integers?

Answer 1

The formula $f$ is the probability mass function for the Poisson distribution. That formula, as explained in the video, can be used to calculate the probability of a given value under the assumed distribution. The related cumulative distribution function $F$ can be used to generate random numbers following the distribution:

1. Use the CDF to partition the interval $(0,1)$ into subintervals: $(0, F(x_1))$, $(F(x_1), F(x_2))$, $etc...$
2. Generate random numbers on the interval $(0,1)$ and see which bin they fall into.

More in this tutorial, which goes through a Poisson example using R. The Poisson PMF and CDF are available in scipy.

Answer 2

The function you link to is a random number generator. It does not return the Poisson distribution, but returns random numbers from a Poisson distribution.

That is, it does exactly what its name suggests - gives you random Poisson variates, not the distribution.

The Poisson probability function is of the form $P(X=x) = \frac{e^{-\lambda} \lambda^x}{x!}\,,\quad 0,1,2,\ldots$,
while the distribution function is $P(X\leq x) =\sum_{i=0}^x \frac{e^{-\lambda} \lambda^i}{i!}\,,\quad 0,1,2,\ldots$.

There are a variety of methods for generating random numbers from this distribution, which will (almost always) begin with a source of uniformly distributed random numbers on $[0,1)$ (notionally continuous, but in practice limited to at best the accuracy with which numbers are represented by the particular implementation on computers).

The scipy function will use one of those methods; which one will be discernable by examining the code (which you'd be better placed to locate than me). However, if I am looking at the right underlying C code that numpy uses (source here), then it uses two different algorithms, depending on the Poisson parameter:

long rk_poisson(rk_state *state, double lam)
{
    if (lam >= 10)
    {
        return rk_poisson_ptrs(state, lam);
    }
    else if (lam == 0)
    {
        return 0;
    }
    else
    {
        return rk_poisson_mult(state, lam);
    }
}

The code for those two functions (rk_poisson_ptrs and rk_poisson_mult) is in the same file, immediately above the quoted code.

Answer 3

I generally use R so my answer here is based on a quick web search. It looks like numpy supports generating random samples from a Poisson distribution and doesn't have functions for computing the probability mass function (PMF) described by the Poisson formula to which you refer. Generating random samples from a distribution can be very useful but as you point out is not the same as computing the PMF which is what you'd need to do to solve the "customer" problem.

It seems like you should be looking at scipy which seems to support the generation of PMF's for a large variety of distributions including Poisson.