#18 Undocumented limitations Poisson implementation

[nlang]

The current Poisson implementation fails for large (>> 100) lambda values, due to numeric intractability of the term 'exp(-lambda)'. Apart from that, the implementation is also slow for large 'lambda', as it runs a 'micro-simulation' which uses about 'lambda' iterations. Technically, this is not a software bug, but a known (but currently undocumented) limitation of the selected implementation.

Suggested quick-fix: document the limits of the selected implementation (they are probably precisely formulated in Shannon, Systems Simulation, 1975, p. 359). For lambda up to 100, I find the results quite accurate (in terms of mean, standard deviation, compared to R output). Users requiring larger lambda values may be suggested to use a normal dist approximation, using mu=s^2=lambda.

Suggested fix: implement a normal approximation for large lambda values, e.g. following 'Ahrens, J.H. and Dieter, U. (1982). Computer generation of Poisson deviates from modified normal distributions. ACM Trans. Math. Software 8, 163-179'. This is for example what is implemented in the R statistic environment.