[R] Re: Generating routine for Poisson random numbers

[Paul Meagher]

> I do think that in an article you should also point out to people that
> there is a lot of numerical code available out there, written by people
> who know a lot more than we do about what they are doing. It's often
> easier than writing your own code and the results are better. One
> advantage of an object-oriented approach is that you can just rip out your
> implementation and slot in a new one if it is better.

Exactly.  That is another reason I do not feel the need to implement the RNG
algorithm perfectly.  If someone really wants a more fool proof rpois-like
algorithm with better running time characteristics they can reimplement
method using the rpois.c normal deviates approach.

BTW, here is what my final Poisson RNG method looks like coded in PHP - it
is modelled after the JSci library approach.  I am only showing the
constructor and the RNG method:

class PoissonDistribution extends ProbabilityDistribution {

  var $lambda;

  function PoissonDistribution($lambda=1) {
    if($lambda <= 0.0) {
      die("Lamda parameter should be positive.");
    }
    $this->lambda = $interval;
  }

  function RNG($num_vals=1) {
    if ($num_vals < 1) {
      die("Number of random numbers to return must be 1 or greater");
    }
    for ($i=0; $i < $num_vals; $i++) {
      $temp  = 0;
      $count = -1;
      while($temp <= $this->lambda) {
        $rand_val = mt_rand() / mt_getrandmax();
        $temp      = $temp - log($rand_val);
        $count++;
      }
      // record count value(s)
      if ($num_vals == 1) {
        $counts = $count;
      } else {
        $counts[$i] = $count;
      }
    }
    return $counts;
  }

My simple eyeball tests indicate that the algorithm appears to generate
unbiased estimates of the expected mean and variance given lambas in the
range of .02 and 900.  I guess to confirm the unbiasedness I would need to
generate a bunch of sample estimates of the mean and variance from my
Poisson random number sequences, plot the relative frequency of these
estimates, and see if the central tendency of the estimates correspond to
the mean and variance expected theoretically for a poisson random variable
(i.e., mean and variance = lambda).

I see why the performance characteristics get bad when lambda is big - the
counting process involves more iterations.  Most of the text book examples
never use a lambda this big, often lambda is less than 100 and often not
less than .02 or so.  In other words, the typical parameter space for the
algorithm may be such that areas where it breaks down are not that common in
practice.

I think this will be a perfectly acceptable RNG for a Poisson random
variable provided you don't use unusually large or small lambda values - if
I knew the break down range, I could implement a check-range test to
disallow usage of the function for that range.

Not sure yet exactly what characteristics of the algorithm would lead it to
behave incorrectly at extremely small or large lambda values?

BTW, is this simple method of generating a poisson random number discussed
in detail in any other books or papers that I might consult?

Regards,
Paul Meagher