[Numpy-discussion] problem with poisson generators

[Bruce Southey]

Hi,
What is the seed used? 
It is really important that you set the seed?
Did you build Python and numarray from source?
Can you reduce your code to a few lines that have the problem? 

Ranlib uses static floats so it may relate to numerical precision (as
well as the random uniform random number generator). Really the only
way is for you to debug the ranlib.c file

Note that R provides a standalone math library in C that includes the
random number generators so you could you those instead. SWIG handles
it rather well.

Regards
Bruce

On 7/13/05, Flavio Coelho <fccoelho at gmail.com> wrote:
> Well,
>   I am definetly glad that someone has also stumbled onto the same problem. 
>  
>  But it is nevertheless intriguing, that you can run poisson a million times
> with mean  zero or negative(it assumes zero mean inthis case) without any
> problems by itself. But when I call it within my code, the rate of error is
> very high (I would say it returns a wrong result every time, but I haven't
> checked every answer).
>  
>  Meanwhile, my solution will be:
>  
>  import rpy
>  
>  n = rpy.r.rpois(n,mean)
>  
>  I don't feel I can trust poisson  while this "funny" behavior is still
> there...
>  If someone has any Idea of how I could trace this bug please tell me, and
> I'll hunt it down. At least I can reproduce it in a very specific context. 
>  
>  thanks,
>  
>  Flávio
> 
> 2005/7/12, Sebastian Haase <haase at msg.ucsf.edu>:
> > Hi Flavio!
> > I had reported this long time ago and this list (about numarray).
> > Somehow this got more or less dismissed. If I recall correctly the
> argument
> > was that nobody could reproduce it (I ran this on Debian Woody , py2.2,
> (with
> > CVS numarray at the time).
> > 
> > I ended up writting my own wrapper(s):
> > def poissonArr(shape=defshape, mean=1):
> >         from numarray import random_array as ra
> >         if mean == 0:
> >                 return zeroArrF(shape) 
> >         elif mean < 0:
> >                 raise "poisson not defined for mean < 0"
> >         else:
> >                 return ra.poisson(mean, shape).astype(na.UInt16)
> > 
> > def poissonize(arr):
> >         from numarray import random_array as ra
> >         return na.where(arr<=0, 0, ra.poisson(arr)).astype(na.UInt16)
> > 
> > (I use the astype( na.UInt16) because of some OpenGL code)
> > 
> > Just last week had this problem on a windows98 computer (python2.4).
> > 
> > This should get sorted out ...
> > 
> > Thanks for reporting this problem.
> > Sebastian Haase
> > 
> > 
> > 
> > 
> > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
> > > Hi,
> > >
> > > I am having problems with the poisson random number generators of both
> > > Numarray and Numeric.
> > > I can't replicate it when calling the function from the python cosonle,
> but 
> > > it has consistently malfunctioned when used within one of my scripts.
> > >
> > > What happens is that it frequently return a value greater than zero when
> > > called with zero mean: poisson(0.0)
> > >
> > > Unfortunately My program is too big to send attached but I have
> confirmed
> > > the malfunction by printing both the mean and the result whenever it
> spits
> > > out a wrong result.
> > >
> > > This is very weird indeed, I have run poisson millions of times by itsel
> on 
> > > the python console, without any problems...
> > >
> > > I hope it is some stupid mistake, but when I replace the poisson
> function
> > > call within my program by the R equivalent command (rpois) via the rpy 
> > > wrapper, everything works just fine...
> > >
> > > any Ideas?
> > 
> 
> 
> 
> -- 
> 
> Flávio Codeço Coelho
> registered Linux user # 386432
> ---------------------------
> Great spirits have always found violent opposition from mediocrities. The 
> latter cannot understand it when a man does not thoughtlessly submit to
> hereditary prejudices but honestly and courageously uses his intelligence.
>         Albert Einstein 
>