Re: [Bug-gsl] two bugs gsl_ran_multinomial on degenerate input, suggestions, and working code

[Ralph Silva]

Hi Alan,

Just try Rmath library (a standalone one) from www.r-project.org. I prefer
to use only one library like GSL, but sometimes I use Rmath for random
distributions for comparisons.  It's open source GPLed.

Best,
Ralph.


On 8/31/07, address@hidden <address@hidden> wrote:
>
> Greetings,
>
> It appears gsl_ran_multinomial was coded not to crap out due to incidental
> negatives in the probability dist. due to floating point artifacts, but
> the handling of these terms was inconsistent and could lead to silently
> absurd results.
>
> I didn't at first intend to go this far.
>
> I'm coding an algorithmic simulation in C++ which makes heavy use of
> multinomials and I discovered the other day that the Boost library has
> only a minimal binomial implementation (Bernoulli trials), and that lead
> me to explore the much improved algorithm in GSL.
>
> Perhaps my problem is unusual, but many vertices in my graph have large K
> and small N.  For efficiency, I was hoping to see GSL multinomial
> implement early-out once all the elements are placed, but it does not.
>
> The multinomial algorithm is quite simple, and on my way to a C++
> implementation that suits my own needs, I decided to code up some fixes
> and a proof-of-concept in C, which the GSL project could potentially adopt
> as it sees fit.
>
> I have supplying a working implementation with two bugs fixes, early-out
> working properly, and a small set of unit tests to catch any glaring
> errors.   I compile it under g++ with the command line:
>
> g++ -l gsl -l gslcblas -O3 -Wall gsl_multinomial_patch.c
>
> You are welcome to incorporate my changes into the GSL project under the
> GPL in any way you see fit.  However, I consider my contributions to this
> C version derivative of the C++ version I have not yet finished, which
> will be released under a BSD license, or potentially a Boost license if it
> goes that far.
>
> What follows is output from my test harness for both my own version, and
> the GSL version.  I think I was linking to 1.8, but reading sources for
> 1.9.
>
> Output from my version is prefixed with ++, the GSL version with --.   The
> existing GSL anomolies can be observed in the tests entitled "bad
> addition" and "eschew negatives".
>
> My version returns an index as part of the early-out optimization.  This
> is the length of the non-zero prefix of the returned value.  My test
> harness is not instrumented to count calls to gsl_ran_binomial() or
> determine execution time, but I am quite certain my routine does call
> gsl_ran_binomial() significantly less than the GSL routine for the test
> case demonstrated.
>
> My version is also more careful to establish strong and formal
> post-conditions, which my test harness double checks.
>
> The source code which produced this listing is attached.  My
> gsl_rand_multinomial revision is written in a pure C idiom (I believe),
> but the rest slides into a C-like style with unwitting C++ influences.
>
> If you don't wish to add a return value to the function signature, I
> suggest a future addition of gsl_multinomial_sparse() and hoisting sum
> p[k] to an initialization routine to further alleviate the O(K)
> performance term.
>
> Feel free to contact me with any questions.
>
> Allan Stokes
>
>
> ===Test harness output===
>
> Begin <empty dist: K=(), N=40>
> ++  0:
> --   :
> End <empty dist>
> Begin <null dist: K=(0.000000, 0.000000, 0.000000), N=40>
> ++  0: 0 0 0
> --   : 0 0 0
> End <null dist>
> Begin <bad addition: K=(0.500000, 0.500000, -0.250000), N=100>
> ++  2: 47 53 0
> ++  2: 56 44 0
> ++  2: 43 57 0
> ++  2: 54 46 0
> ++  2: 46 54 0
> ++  2: 51 49 0
> ++  2: 46 54 0
> ++  2: 46 54 0
> ++  2: 48 52 0
> ++  2: 55 45 0
> --   : 69 31 0
> --   : 73 27 0
> --   : 64 36 0
> --   : 68 32 0
> --   : 67 33 0
> --   : 61 39 0
> --   : 68 32 0
> --   : 66 34 0
> --   : 61 39 0
> --   : 63 37 0
> End <bad addition>
> Begin <eschew negatives: K=(0.010000, -100.000000, 0.010000), N=100>
> ++  3: 55 0 45
> ++  3: 53 0 47
> ++  3: 51 0 49
> ++  3: 57 0 43
> ++  3: 52 0 48
> ++  3: 51 0 49
> ++  3: 59 0 41
> ++  3: 43 0 57
> ++  3: 48 0 52
> ++  3: 45 0 55
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> --   : 0 0 100
> End <eschew negatives>
> Begin <early out: K=(0.800000, 0.160000, 0.040000), N=10>
> ++  3: 8 1 1
> ++  3: 6 1 3
> ++  2: 9 1 0
> ++  3: 8 1 1
> ++  1: 10 0 0
> ++  3: 7 2 1
> ++  3: 9 0 1
> ++  2: 9 1 0
> ++  2: 7 3 0
> ++  2: 9 1 0
> ++  3: 8 1 1
> ++  3: 6 3 1
> ++  3: 8 1 1
> ++  3: 7 1 2
> ++  3: 9 0 1
> ++  2: 9 1 0
> ++  3: 6 2 2
> ++  2: 8 2 0
> ++  2: 9 1 0
> ++  2: 9 1 0
> ++  2: 9 1 0
> ++  2: 8 2 0
> ++  2: 8 2 0
> ++  2: 7 3 0
> ++  2: 9 1 0
> ++  3: 6 3 1
> ++  2: 9 1 0
> ++  2: 9 1 0
> ++  2: 9 1 0
> ++  2: 8 2 0
> ++  2: 9 1 0
> ++  3: 7 1 2
> ++  3: 8 1 1
> ++  3: 8 1 1
> ++  3: 8 1 1
> ++  2: 6 4 0
> ++  3: 9 0 1
> ++  3: 8 1 1
> ++  1: 10 0 0
> ++  3: 8 1 1
> --   : 8 1 1
> --   : 8 2 0
> --   : 8 2 0
> --   : 9 1 0
> --   : 8 2 0
> --   : 10 0 0
> --   : 6 3 1
> --   : 8 1 1
> --   : 9 1 0
> --   : 9 1 0
> --   : 8 1 1
> --   : 10 0 0
> --   : 8 2 0
> --   : 8 2 0
> --   : 8 2 0
> --   : 8 2 0
> --   : 8 2 0
> --   : 8 2 0
> --   : 10 0 0
> --   : 9 0 1
> --   : 8 2 0
> --   : 8 2 0
> --   : 8 1 1
> --   : 8 2 0
> --   : 10 0 0
> --   : 7 3 0
> --   : 7 0 3
> --   : 8 1 1
> --   : 6 4 0
> --   : 7 3 0
> --   : 10 0 0
> --   : 7 2 1
> --   : 10 0 0
> --   : 9 1 0
> --   : 8 2 0
> --   : 7 3 0
> --   : 9 1 0
> --   : 7 3 0
> --   : 9 0 1
> --   : 6 2 2
> End <early out>
>
>
>
> _______________________________________________
> Bug-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/bug-gsl
>
>
>

Rmath.h
Description: Binary data