Re: [Numpy-discussion] problem with poisson generators

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

2005/7/13, Bruce Southey <bso...@gm...>:
>=20
> Hi,
> What is the seed used?


I am not setting the seed.

It is really important that you set the seed?


No.

Did you build Python and numarray from source?



Yes, I use Gentoo. I build everything from source.=20

Can you reduce your code to a few lines that have the problem?



Not really, poisson and binomial are the only two functions from Numeric=20
that I use but they are called inside a rather complex object oriented code=
=20
environment (objects within objetcs, being called recursively...) Bellow is=
=20
the function within which poisson is called:

def stepSEIR_s(self,inits=3D(0,0,0),par=3D(0.001,1,1,0.5,1,0,0,0),theta=3D0=
,=20
npass=3D0,dist=3D'poisson'):
 """
 Defines an stochastic model SEIR:
 - inits =3D (E,I,S)
 - par =3D (Beta, alpha, E,r,delta,B,w,p) see docs.
 - theta =3D infectious individuals from neighbor sites
 """
 E,I,S =3D inits
 N =3D self.parentSite.totpop
 beta,alpha,e,r,delta,B,w,p =3D par
 Lpos_esp =3D beta*S*((I+theta)/(N+npass))**alpha #Number of new cases
=20
 if dist =3D=3D 'poisson':
 Lpos =3D poisson(Lpos_esp)
## if theta =3D=3D 0 and Lpos_esp =3D=3D 0 and Lpos > 0:=20
## print Lpos,Lpos_esp,S,I,theta,N,self.parentSite.sitename
 elif dist =3D=3D 'negbin':
 prob =3D I/(I+Lpos_esp) #convertin between parameterizations
 Lpos =3D negative_binomial(I,prob)
 self.parentSite.totalcases +=3D Lpos #update number of cases=20
 Epos =3D (1-e)*E + Lpos
 Ipos =3D e*E + (1-r)*I
 Spos =3D S + B - Lpos=20
 Rpos =3D N-(Spos+Epos+Ipos)
 #self.parentSite.totpop =3D Spos+Epos+Ipos+Rpos
 self.parentSite.incidence.append(Lpos)
 if not self.parentSite.infected:
 if Lpos > 0:
 self.parentSite.infected =3D self.parentSite.parentGraph.simstep
 self.parentSite.parentGraph.epipath.append((
self.parentSite.parentGraph.simstep,self.parentSite,self.parentSite.infecto=
r
))

=20
 return [Epos,Ipos,Spos]


I wonder if called by itself it would trigger the problem... The commented=
=20
Lines Is where I catch the errors: when poisson returns a greater than zero=
=20
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.=20

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?=20

thanks Bruce,

Fl=E1vio

Regards
> Bruce
>=20
> On 7/13/05, Flavio Coelho <fcc...@gm...> wrote:
> > Well,
> > I am definetly glad that someone has also stumbled onto the same=20
> problem.
> >
> > But it is nevertheless intriguing, that you can run poisson a million=
=20
> times
> > with mean zero or negative(it assumes zero mean inthis case) without an=
y
> > problems by itself. But when I call it within my code, the rate of erro=
r=20
> is
> > very high (I would say it returns a wrong result every time, but I=20
> 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,=
=20
> and
> > I'll hunt it down. At least I can reproduce it in a very specific=20
> 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=
=20
> both
> > > > Numarray and Numeric.
> > > > I can't replicate it when calling the function from the python=20
> cosonle,
> > but
> > > > it has consistently malfunctioned when used within one of my=20
> scripts.
> > > >
> > > > What happens is that it frequently return a value greater than zero=
=20
> 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 i=
t
> > spits
> > > > out a wrong result.
> > > >
> > > > This is very weird indeed, I have run poisson millions of times by=
=20
> 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=
=20
> 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.=
=20
> The
> > latter cannot understand it when a man does not thoughtlessly submit to
> > hereditary prejudices but honestly and courageously uses his=20
> intelligence.
> > Albert Einstein
> >
>=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