Mailman 3 RE: [Numpy-discussion] problem with poisson generators - NumPy-Discussion

July 15, 2005

      Hi Flavio,

Could you give us every call to poisson or negative_binomial (ie all functions related to random numbers) preceeded with the seed ?

Adding before your declaration of stepSEIR_s code like:

randoutput = file("randoutput.py", "w")
old_poisson = poisson
def poisson(m):
    print >> randoutput, "print get_seed(), poisson(%s)"%m
    result = old_poisson(m)
    print >> randoutput, "# result = %s"%result
    return result

old_negative_binomial = negative_binomial
def negative_binomial(i,p):
    print >> randoutput, "print get_seed(), negative_binomial(%s,%s)"%(i,p)
    result = old_negative_binomial(i,p)
    print >> randoutput, "# result = %s"%result
    return result

will ouput every call.
We can then easily try on our machines what it gives.

Best,

Sebastien
...
-----Original Message-----
From: numpy-discussion-admin@lists.sourceforge.net
[mailto:numpy-discussion-admin@lists.sourceforge.net]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
On Wed, 2005-07-13 at 12:13 -0300, Flavio Coelho wrote:
...
2005/7/13, Bruce Southey <bsouthey@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
...
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
complex object
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.par
entSite.infector))
...
return [Epos,Ipos,Spos]
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@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@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
-------------------------------------------------------
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=7477&alloc_id=16492&op=click
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/numpy-discussion
=======================================================
This message is confidential. It may also be privileged or otherwise protected 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 sent to and from Electrabel may be monitored to ensure compliance with internal policies and to protect our business. Emails are not secure and cannot be 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
=======================================================

RE: [Numpy-discussion] problem with poisson generators

Sebastien.deMentendeHorne＠electrabel.com

tags

participants (1)