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

[allan]

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>

gsl_multinomial_patch.c
Description: Text Data