[Numpy-discussion] problem with poisson generators

Todd Miller jmiller at stsci.edu
Fri Jul 15 18:52:31 CDT 2005


On Fri, 2005-07-15 at 21:30 -0300, Flavio Coelho wrote:
> 
> 
> 2005/7/15, Todd Miller <jmiller at stsci.edu>:
> 
> Todd,
> 
> the  function  within which the problem happens runs fine by itself,
> as well as the simplified version you sent me.

Maybe you could stick in an assert rather than the print and then use
pdb.pm() to find the actual parameters leading to the failure.  (The
parameters may or may not matter, i.e. the problem may be hidden state
somewhere your program accumulates in some complicated way,  but maybe
we'll get lucky...)

> I am running the code on a AMD64 (though all my software and OS is
> compiled to 686 - Gentoo Linux). Is there any issue I should be aware
> of,  regarding this specific cpu? 

Other than the fact that 64-bit addressing functionality isn't where we
want it to be,  there are no remaining issues I'm aware of for AMD64.
But the addressing is currently hamstrung by 32-bit limits in Python's
sequence protocol so you should be aware that even with your AMD64 you
may run out of head room.

> I can try to run the software on another machine...

I wouldn't bother:  our goal should be to replicate the problem simply
and that suggests to me maintaining identical initial conditions.  BTW,
I'm using AMD64 too so that won't be another variable.

Regards,
Todd

> 
>         Having this code is a good start to solving the problem.  I
>         think the
>         next step is to simplify your example to make it runnable and
>         provide 
>         known inputs for all the parameters which lead to a failure
>         for you.
>         
>         Being really literal (spelling out every single damn thing)
>         cuts down on
>         speculation and inefficiency on our end.
>         
>         It would also be good to express the condition (or it's
>         inverse) you 
>         think is an error as an assertion,  so something like this
>         might be what
>         we need:
>         
>         from numarray.random_array import poisson
>         
>         E, I, S = (0,0,0)
>         beta,alpha,e,r,delta,B,w,p = (0.001,1,1,0.5,1,0,0,0)
>         theta = 0
>         npass = 0
>         N = 100 # total guess here for me
>         Lpos_esp = beta*S*((I+theta)/(N+npass))**alpha #Number of new
>         cases
>         Lpos = poisson(Lpos_esp)
>         assert not (theta == 0 and Lpos_esp == 0 and Lpos > 0)
>         
>         The problem is,  the above works for me.  Make it sane and get
>         it to
>         expose the error for you and I'll look into this some more.
>         
>         Regards,
>         Todd
>         
>         >
>         > I wonder if called by itself  it would trigger the
>         problem... The 
>         > commented Lines Is where I catch the errors: when poisson
>         returns a
>         > greater than zero number, when called with mean 0.
>         >
>         >
>         >
>         >         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
>         >
>         > I'll see what I can do. 
>         >
>         >
>         >         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.
>         >
>         >
>         > I think thats what Rpy already does...is it not?
>         >
>         > thanks Bruce,
>         >
>         > Flávio
>         >
>         >
>         >         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
>         >         >
>         >
>         >
>         >
>         > --
>         > 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
>         
> 
> 
> 
> -- 
> 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 





More information about the Numpy-discussion mailing list