pioneerspacesim/pioneer — pull request #3242 by impaktor (2014-11-07 07:52:48+00:00)
Title: 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