0

I'm trying to generate multinomial random variables as fast as possible. And I learned that gsl_ran_multinomial could be a good choice. However, I tried to use it based on the answer in this post: https://stackoverflow.com/a/23100665/21039115, and the results were always wrong.

In detail, my code is

// [[Rcpp::export]]
arma::ivec rmvn_gsl(int K, arma::vec prob) {
    gsl_rng *s = gsl_rng_alloc(gsl_rng_mt19937); // Create RNG seed
    arma::ivec temp(K);
    gsl_ran_multinomial(s, K, 1, prob.begin(), (unsigned int *) temp.begin());
    gsl_rng_free(s); // Free memory
    return temp;
}

And the result was something like

rmvn_gsl(3, c(0.2, 0.7, 0.1))
     [,1]
[1,]    1
[2,]    0
[3,]    0

which is ridiculous.

I was wondering if there was any problem exist in the code... I couldn't find any other examples to compare. I appreciate any help!!!

UPDATE:

I found the primary problem here is that I didn't set a random seed, and it seems that gsl has its own default seed (FYI: https://stackoverflow.com/a/32939816/21039115). Once I set the seed by time, the code worked. But I will go with rmultinom since it can even be faster than gsl_ran_multinomial based on microbenchmark.

Anyway, @Dirk Eddelbuettel provided a great example of implementing gsl_ran_multinomial below. Just pay attention to the random seeds issue if someone met the same problem as me.

dywu0106
  • 3
  • 2
  • Welcome to StackOverflow. It's a little hard to start on your question as there are multiple things wrong all at once. As the saying goes, "learn to walk before you try to run". Maybe you should focus on the other question and _try to just create three random draws_ and maybe print them. Then maybe try to return them. Mixing `arma` and `gsl` vectors is not a good idea, and passing double vectors cast to `int` is a terrible idea. Decompose. check indiividual subtasks. Good luck! – Dirk Eddelbuettel Jan 18 '23 at 21:20
  • Thank you so much!!! You are right, I wrongly passed double vectors cast to int, and I've adjusted the code. But I still want to try to keep both `arma` and `gsl`. `arma` is a great part of my other functions, and `gsl` can provide an efficient random variable generation. – dywu0106 Jan 18 '23 at 22:14
  • There is no reason not to _eventually_ but the GSL framework (in C only) does not know or work with Armadillo. So I would use RcppGSL to get yourself a vector, and then make an arma vector from it. – Dirk Eddelbuettel Jan 18 '23 at 22:15
  • Also see our `Rmath.h` header -- we expose the C level function of R in `inline void rmultinom(int n, double* prob, int k, int* rn)`. Might be easier to deal with.... – Dirk Eddelbuettel Jan 18 '23 at 22:49

1 Answers1

1

Here is a complete example taking a double vector of probabilities and returning an unsigned integer vector (at the compiled level) that is mapped to an integer vector by the time we are back in R:

Code

#include <RcppGSL.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

// [[Rcpp::depends(RcppGSL)]]

// [[Rcpp::export]]
std::vector<unsigned int> foo(int n, const std::vector <double> p) {
    int k = p.size();
    std::vector<unsigned int> nv(k);
    gsl_rng_env_setup();
    gsl_rng *s = gsl_rng_alloc(gsl_rng_mt19937); // Create RNG instance
    gsl_ran_multinomial(s, k, n, &(p[0]), &(nv[0]));
    gsl_rng_free(s);

    return nv;
}


/*** R
foo(400, c(0.1, 0.2, 0.3, 0.4))
*/

Output

> Rcpp::sourceCpp("~/git/stackoverflow/75165241/answer.cpp")

> foo(400, c(0.1, 0.2, 0.3, 0.4))
[1]  37  80 138 145
> 
Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725
  • Hi Dirk, I really appreciate it! I would just like to set a total number `n = 1`. If I did that and repeated it several times, the results wouldn't change and always be `c(1, 0, 0, 0)` even with the probabilities `c(0.1, 0.2, 0.3, 0.4)`. Then I realized that was because of the default seed `gsl` used. Once I changed the seeds by time, everything worked well. But I have to say `rmultinom` is much easier, and it seems that it can even be faster than `gsl_ran_multinomial` based on microbenchmark! Thanks again for your kind help!!! :) – dywu0106 Jan 19 '23 at 01:07