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, [EMAIL PROTECTED] <[EMAIL PROTECTED]> 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 > [email protected] > http://lists.gnu.org/mailman/listinfo/bug-gsl > > >
Rmath.h
Description: Binary data
_______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
