[Numpy-discussion] problem with poisson generators

Wed Jul 13 08:17:21 EDT 2005

2005/7/13, Bruce Southey <bsouthey at gmail.com>:
> 
> 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. 

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 
that I use but they are called inside a rather complex object oriented code 
environment (objects within objetcs, being called recursively...) Bellow is 
the function within which poisson is called:

def stepSEIR_s(self,inits=(0,0,0),par=(0.001,1,1,0.5,1,0,0,0),theta=0, 
npass=0,dist='poisson'):
 """
 Defines an stochastic model SEIR:
 - inits = (E,I,S)
 - par = (Beta, alpha, E,r,delta,B,w,p) see docs.
 - theta = infectious individuals from neighbor sites
 """
 E,I,S = inits
 N = self.parentSite.totpop
 beta,alpha,e,r,delta,B,w,p = par
 Lpos_esp = beta*S*((I+theta)/(N+npass))**alpha #Number of new cases

 if dist == 'poisson':
 Lpos = poisson(Lpos_esp)
## if theta == 0 and Lpos_esp == 0 and Lpos > 0: 
## print Lpos,Lpos_esp,S,I,theta,N,self.parentSite.sitename
 elif dist == 'negbin':
 prob = I/(I+Lpos_esp) #convertin between parameterizations
 Lpos = negative_binomial(I,prob)
 self.parentSite.totalcases += Lpos #update number of cases 
 Epos = (1-e)*E + Lpos
 Ipos = e*E + (1-r)*I
 Spos = S + B - Lpos 
 Rpos = N-(Spos+Epos+Ipos)
 #self.parentSite.totpop = Spos+Epos+Ipos+Rpos
 self.parentSite.incidence.append(Lpos)
 if not self.parentSite.infected:
 if Lpos > 0:
 self.parentSite.infected = self.parentSite.parentGraph.simstep
 self.parentSite.parentGraph.epipath.append((
self.parentSite.parentGraph.simstep,self.parentSite,self.parentSite.infector
))

 return [Epos,Ipos,Spos]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20050713/a3c88678/attachment.html>