Re: [Numpy-discussion] problem with poisson generators

[Flavio C. <fcc...@gm...>]

2005/7/15, Todd Miller <jm...@st...>:
Todd,

the function within which the problem happens runs fine by itself, as well=
=20
as the simplified version you sent me.

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

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.
>=20
> Being really literal (spelling out every single damn thing) cuts down on
> speculation and inefficiency on our end.
>=20
> 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:
>=20
> from numarray.random_array import poisson
>=20
> E, I, S =3D (0,0,0)
> beta,alpha,e,r,delta,B,w,p =3D (0.001,1,1,0.5,1,0,0,0)
> theta =3D 0
> npass =3D 0
> N =3D 100 # total guess here for me
> Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha #Number of new cases
> Lpos =3D poisson(Lpos_esp)
> assert not (theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0)
>=20
> 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.
>=20
> Regards,
> Todd
>=20
> >
> > 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=E1vio
> >
> >
> > Regards
> > Bruce
> >
> > On 7/13/05, Flavio Coelho <fcc...@gm...> 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 =3D 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=E1vio
> > >
> > > 2005/7/12, Sebastian Haase <ha...@ms...>:
> > > > 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=3Ddefshape, mean=3D1):
> > > > from numarray import random_array as ra
> > > > if mean =3D=3D 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<=3D0, 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=E1vio Code=E7o 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=E1vio Code=E7o 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
>=20
>=20


--=20
Fl=E1vio Code=E7o 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