libmir / mir-random

Advanced Random Number Generators
http://mir-random.libmir.org/
32 stars 15 forks source link

Implementation of Multinomial random variable #118

Closed fccoelho closed 4 years ago

fccoelho commented 4 years ago

I have ported this from GSL as described on issue #117. it has been adapted to Mir-random API and has a simple unittest which are passing.

There is a minor deviation from the way other ndvariables in that the opCall is not void, but rather returns the result expected. I try but could not make it alter the result variable passed in in the outer scope, even though my implementation follows exactly the other variable's implementations. The only difference is that multinomial returns an array of long integers instead of doubles. If there's a way to fix this, please let me know and I'll change the code, but as it is the function works correctly. It is worth pointing out that the opCall of the scalar random variables, defined in variable.d, are also not void, returning the results directly.

9il commented 4 years ago

GSL likely isn't compatible with BSL-1.0. So, the direct port can't be accepted. However, if you understand the math and implement the algorithm from scratch, then it can be accepted.

fccoelho commented 4 years ago

GSL likely isn't compatible with BSL-1.0. So, the direct port can't be accepted. However, if you understand the math and implement the algorithm from scratch, then it can be accepted.

Regarding this licensing issue, the GSL implementation comes straight from this reference:

 C.S. David, The computer generation of multinomial random variates,
   Comp. Stat. Data Anal. 16 (1993) 205-217

The algorithm is so simple that there is no other way I see to implement it efficiently. I can bet that all other implementations of the RNG are based on the same algorithm.

9il commented 4 years ago

few nitpicks, other LGTM

fccoelho commented 4 years ago

Is there anything left to fix here, before it can be merged?

9il commented 4 years ago

norm member should be removed as we assert in the constructor the norm is about 1., opCall shouls just use 1 as well.