[SciPy-User] NumPy Binomial BTPE method Problem

[Josh Lawrence]

Hey all,

I received access to the paper and it seems it was originally based
purely on the paper written by Kachitvichyanukul in 1988. I still
think there's a whoopsies with the if ... else if ... else, block
though.

On Wed, Oct 3, 2012 at 11:54 AM, Josh Lawrence
<josh.k.lawrence at gmail.com> wrote:
> 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.



-- 
Josh Lawrence