Random number from Poisson distribution

# Random number from Poisson distribution

## Why?

Whenever we have a series of independent random events of integer values (like radioactive decay in some time interval), where we know their average, their distribution follows a Poisson distribution, plus I think I need it. For instance, this can be used to determine how many random events have taken place, rather than simulating each separately by drawing a uniform random number for each.

## It works

Demonstration of code output. Poisson distribution only takes the average (usuallt called `lambda`) as input parameter:

![poisson_3](https://cloud.githubusercontent.com/assets/619390/4949907/276d1e9a-6652-11e4-9519-1b24e71ade4e.png)

![poisson_15](https://cloud.githubusercontent.com/assets/619390/4949913/35429234-6652-11e4-937e-d08ebe898a04.png)


## Notes on implementation

There are quite a number of ways to implement this. The method used here is from Knuth, I think. It is a lot faster than the naive cumulative probability way. There are faster methods still, which scale better for high lambda, but we might as well use the `Random::Normal` for that. Those algorithms are of interest when one needs a huge amount of random Poisson distributed numbers. If that is to happen in the future, I think [this](http://www.evolvedmicrobe.com/Blog.html) ([alternative link](http://evolvedmicrobe.com/blogs/?p=6)) will be a gold mine of the different methods out there.

---

## Appendix:

### Gnuplot script used above

Because last time, #2651, I was asked to provide my plotting script I think, through IRC. Plus I now had to rewrite it.

    set terminal png enhanced
    set xlabel "Random number"
    set ylabel "freq"
    set title "Poisson distribution"
    set yrange[0:]

    bin(x,width)=width*floor(x/width)
    binwidth = .1

    f(k,lambda) = lambda**k * exp(-lambda) / gamma(k+1)
    set output "poisson_3.png"
    plot "~/3.dat" using (bin($1,binwidth)):(1.0) smooth freq with boxes ti "{/Symbol l} = 3",\
         f(x,3)*5000 ti "{/Symbol l} = 3"

    set output "poisson_15.png"
    plot "~/15.dat" using (bin($1,binwidth)):(1.0) smooth freq with boxes ti "{/Symbol l} = 15",\
         f(x,15)*5000 ti "{/Symbol l} = 15"


### Generating testing data

I put this in a file `data/modules/foo.lua` to generate data used above.

    local Engine = import("Engine")
    for i=1,5000 do
    	print(Engine.rand:Poisson(3))
    end