[R] Poisson generation (was Viewing function source)

[Prof Brian Ripley]

See any good book on simulation: one from 1987 springs to mind.

The hard issue is to find a method that works near optimally for all
values of lambda, especially moderately large ones.  No simple algorithm
does that, and for something like R, the issue is not compactness but good 
worst-case behaviour.

I do think you will learn more by reading a book than examples of code, at 
least until you know the monographs backwards (and forwards).

On Tue, 26 Aug 2003, Paul Meagher wrote:

> > Just type the name of the function to see the R code
> >   > rpois
> >   function (n, lambda)
> >   .Internal(rpois(n, lambda))
> >
> > But in this case it tells you that rpois is implemented in C code :(
> >
> >  By convention, it is likely to be a function called do_rpois,
> >  however in this case we have an exception to the convention.  You can
> > look in src/main.names.c for the name of the C function.
> >
> > abacus% fgrep rpois names.c
> > {"rpois",       do_random1,     3,      11,     2,      {PP_FUNCALL,
> > PREC_FN,0}},
> >
> > so the function do_random1 is actually involved. If you grep for that, you
> > are directed to src/main/random.c and you will eventually find that the
> > underlying code is in src/nmath/rpois.c and is commented with a reference
> > to the ACM Transactions on Mathematical Software.
> 
> Thank you for helping me with the research on this.
> 
> I am interested in the rpois source in particular because I am wanting to
> implement a compact algorithm for it.
> 
> The rpois.c code is estimated at around 140 lines of code when you subtract
> out line spaces and commenting.   Fairly inscrutable if you haven't read the
> article or don't have a mathematical programming background.
> 
> I had hoped that I would see a compact algorithm of around 20 lines of code
> based on a counting process sensitive to a lambda parameter.  Instead it
> uses a normal deviates method that is probably much more optimal in many
> mathematical senses.
> 
> Personally, I would like to see a counting process implementation of an
> poisson random number generator.  I suspect it would be much slower than
> rpois.c (because it would likely depend upon setting a num_frames iteration
> counter) and less accurate, but would be much more compact and give more
> intuitive insight and understanding of a numerical process that generates

One of the possible processes.

> the poisson random deviates.  I am not very fluent yet in R programming to
> attempt this.  I could take a stab at it with PHP which I am much more
> fluent in if I am not barking up the wrong tree on this conjecture.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595