testing randp

[Michael Creel]

Hi, that last test I sent in was using poisson_rnd - sorry I didn't read 
carefully enough to see that there was a new generator. The following 
function can be used to test that it's working ok. Since the Poisson mean and 
variance are both lambda, evidence that the sample moment divided by lambda 
is different from 1 would cast doubt upon the generator. I find that it seems 
to work fine for all lambdas except very big ones, where it segfaults (1e8 is 
ok, 1e9 will crash it).

function test_randp(lambda, reps)
 tic();
 y = randp(lambda, reps, 1);
 printf("lambda %f\n",lambda)
 printf("time to generate %d random draws: %f\n", reps, toc());
 printf("sample mean / true mean %f\n", mean(y)/lambda);
 printf("sample variance / true variance %f\n", var(y)/lambda);
endfunction

The test could be formalized (just using a t-test that the sample moment 
equals lambda), and higher moments could be used, but I ain't got the energy.

Michael

> on 10/6/04 8:53 PM, Paul Kienzle at address@hidden wrote:
> > I've posted randp() to octave-forge/FIXES.   The code for
> > the 5 (?!) generators  is in randpoisson.c.  There are
> > different generators for large (>10) and small lambda (<10),
> > and for repeated and varying lambda.  Some algorithms have
> > considerable initialization overhead when changing
> > lambda.   Also, for huge lambda (>1e8) a normal approximation
> > is used.
> >
> > Be careful to test that the generated numbers look
> > reasonable for all these cases before putting your
> > trust in this code.  Test showing that the resulting
> > numbers are indeed Poisson would be appreciated.



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------