Menu
▾
▴
RE: [Numpy-discussion] problem with poisson generators
|
From: <Seb...@el...> - 2005-07-15 11:07:56
|
Hi Flavio,
Could you give us every call to poisson or negative_binomial (ie all functi=
ons related to random numbers) preceeded with the seed ?
Adding before your declaration of stepSEIR_s code like:
randoutput =3D file("randoutput.py", "w")
old_poisson =3D poisson
def poisson(m):
print >> randoutput, "print get_seed(), poisson(%s)"%m
result =3D old_poisson(m)
print >> randoutput, "# result =3D %s"%result
return result
old_negative_binomial =3D negative_binomial
def negative_binomial(i,p):
print >> randoutput, "print get_seed(), negative_binomial(%s,%s)"%(i,p)
result =3D old_negative_binomial(i,p)
print >> randoutput, "# result =3D %s"%result
return result
will ouput every call.
We can then easily try on our machines what it gives.
Best,
Sebastien
> -----Original Message-----
> From: num...@li...
> [mailto:num...@li...]On Behalf Of Todd
> Miller
> Sent: vendredi 15 juillet 2005 12:54
> To: Flavio Coelho
> Cc: Bruce Southey; Sebastian Haase; numpy-discussion
> Subject: Re: [Numpy-discussion] problem with poisson generators
>=20
>=20
> On Wed, 2005-07-13 at 12:13 -0300, Flavio Coelho wrote:
> >=20
> >=20
> > 2005/7/13, Bruce Southey <bso...@gm...>:
> > Hi,
> > What is the seed used?
> >=20
> > I am not setting the seed.
> > =20
> >=20
> > It is really important that you set the seed?
> >=20
> > No.
> >=20
> >=20
> > Did you build Python and numarray from source?
> >=20
> >=20
> > Yes, I use Gentoo. I build everything from source.=20
> >=20
> >=20
> > Can you reduce your code to a few lines that have=20
> the problem?
> >=20
> >=20
> > Not really, poisson and binomial are the only two functions from
> > Numeric that I use but they are called inside a rather=20
> complex object
> > oriented code environment (objects within objetcs, being called
> > recursively...) Bellow is the function within which poisson=20
> is called:
> >=20
> > def=20
> stepSEIR_s(self,inits=3D(0,0,0),par=3D(0.001,1,1,0.5,1,0,0,0),theta=3D0,
> > 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=20
> #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=20
> parameterizations
> > Lpos =3D negative_binomial(I,prob)
> > self.parentSite.totalcases +=3D Lpos #update number=20
> 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
> >=20
> ((self.parentSite.parentGraph.simstep,self.parentSite,self.par
> entSite.infector))
> >=20
> > =20
> > return [Epos,Ipos,Spos]
> >=20
>=20
> 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
>=20
> Being really literal (spelling out every single damn thing)=20
> 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=20
> 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
> >=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.
> >=20
> >=20
> >=20
> > Ranlib uses static floats so it may relate to numerical
> > precision (as
> > well as the random uniform random number generator). Really
> > the only=20
> > way is for you to debug the ranlib.c file
> >=20
> > I'll see what I can do.=20
> >=20
> >=20
> > Note that R provides a standalone math library in C that
> > includes the
> > random number generators so you could you those=20
> instead. SWIG
> > handles=20
> > it rather well.
> >=20
> >=20
> > I think thats what Rpy already does...is it not?=20
> >=20
> > thanks Bruce,
> >=20
> > Fl=C3=A1vio
> > =20
> >=20
> > 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 problem.
> > >
> > > But it is nevertheless intriguing, that you can=20
> run poisson
> > a million times=20
> > > 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=20
> every time,
> > but I haven't=20
> > > 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=20
> > > 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=C3=A1vio=20
> > >
> > > 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=20
> > > 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):=20
> > > > def poissonArr(shape=3Ddefshape, mean=3D1):
> > > > from numarray import random_array as ra
> > > > if mean =3D=3D 0:
> > > > return zeroArrF(shape)
> > > > elif mean < 0:=20
> > > > 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,=20
> ra.poisson(arr)).astype
> > ( na.UInt16)
> > > >
> > > > (I use the astype( na.UInt16) because of some=20
> OpenGL code)
> > > >
> > > > Just last week had this problem on a windows98 computer
> > (python2.4).
> > > >
> > > > This should get sorted out ...=20
> > > >
> > > > Thanks for reporting this problem.
> > > > Sebastian Haase
> > > >
> > > >
> > > >
> > > >
> > > > On Tuesday 12 July 2005 13:32, Flavio Coelho wrote:
> > > > > Hi,=20
> > > > >
> > > > > I am having problems with the poisson random number
> > generators of both
> > > > > Numarray and Numeric.
> > > > > I can't replicate it when calling the=20
> function from the
> > python cosonle,=20
> > > 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=20
> attached but
> > I have
> > > confirmed
> > > > > the malfunction by printing both the mean and=20
> the result
> > whenever it
> > > spits=20
> > > > > 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=20
> replace the
> > poisson
> > > function
> > > > > call within my program by the R equivalent command
> > (rpois) via the rpy
> > > > > wrapper, everything works just fine...=20
> > > > >
> > > > > any Ideas?
> > > >
> > >
> > >
> > >
> > > --
> > >
> > > Fl=C3=A1vio Code=C3=A7o Coelho
> > > registered Linux user # 386432
> > > ---------------------------
> > > Great spirits have always found violent opposition from
> > mediocrities. The=20
> > > latter cannot understand it when a man does not
> > thoughtlessly submit to
> > > hereditary prejudices but honestly and=20
> courageously uses his
> > intelligence.
> > > Albert Einstein
> > >
> >=20
> >=20
> >=20
> > --=20
> > Fl=C3=A1vio Code=C3=A7o Coelho
> > registered Linux user # 386432
> > ---------------------------
> > Great spirits have always found violent opposition from=20
> mediocrities.
> > The
> > latter cannot understand it when a man does not thoughtlessly submit
> > to=20
> > hereditary prejudices but honestly and courageously uses his
> > intelligence.
> > Albert Einstein=20
>=20
>=20
>=20
> -------------------------------------------------------
> SF.Net email is sponsored by: Discover Easy Linux Migration Strategies
> from IBM. Find simple to follow Roadmaps, straightforward articles,
> informative Webcasts and more! Get everything you need to get up to
> speed, fast. http://ads.osdn.com/?ad_id=3D7477&alloc_id=3D16492&op=3Dclick
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>=20
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D
This message is confidential. It may also be privileged or otherwise protec=
ted by work product immunity or other legal rules. If you have received it =
by mistake please let us know by reply and then delete it from your system;=
you should not copy it or disclose its contents to anyone. All messages se=
nt to and from Electrabel may be monitored to ensure compliance with intern=
al policies and to protect our business. Emails are not secure and cannot b=
e guaranteed to be error free as they can be intercepted, amended, lost or =
destroyed, or contain viruses. Anyone who communicates with us by email is =
taken to accept these risks.
http://www.electrabel.be/homepage/general/disclaimer_EN.asp
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D
|
