[SciPy-User] NumPy Binomial BTPE method Problem

[Josh Lawrence]

Hello all,

I am implementing a binomial random variable in MATLAB. The default
method in the statistics toolbox is extremely slow for large
population/trial size. I am needing to do trials for n as large as
2**28. I found in NumPy some code that implements a binomial random
draw in numpy/random/mtrand/distributions.c. I was trying to convert
the code to MATLAB and the BTPE method seems to have an error in lines
337-341 of distributions.c. The if ... else if ... else statement I
think is incorrect. I think it should be an if ... else ... statement
followed by the contents of the original else which starts on line
337.

The if ... else if ... else block is as follows:

#### begin code snippet ####
if (m < y)
{
    for (i=m; i<=y; i++)
    {
        F *= (a/i - s);
    }
}
else if (m > y)
{
    for (i=y; i<=m; i++)
    {
        F /= (a/i - s);
    }
}
else
{
    if (v > F) goto Step10;
    goto Step60;
}
#### end code snippet ####

>From what I can tell, the variable F is only used in the comparison
within the else{} statment (i.e. the if(v > F) statement) and nowhere
else within the scope of the function.

I also found a fortran implementation here:
http://wstein.org/home/wstein/www/home/mhansen/spkgs_in_progress/octave-3.2.4/src/libcruft/ranlib/ignbin.f
and it appears this is from where the code was originally adapted as
the variable names are the same.

My parsing of fortran GOTOs is a bit rusty, but I think the contents
of the else block in above snippet should be not be conditional.

I don't understand the underlying algorithm very well and don't have
access the the BTPE paper, so I can't comment on the validity of the
fortran code. There just seems to be an error in logic in the above
code. So please have someone who understands it look at it. It appears
Robert Kern wrote the function a decent portion of the file at some
point in the past.

I hope this helps.

Cheers,

--
Josh Lawrence


P.S. I apologize if my email is inconvenient, but I could not figure
out how to tell gmail to set the reply-to field to be
scipy-user at scipy.org.