Using rmultinom with Rcpp

Question

I'd like to use the R function rmultinom in c++ code to be used with Rcpp. I get an error about not enough arguments - I am unfamiliar with what these arguments ought to be, as they do not corresond to the arguments used with the function in R. I also didn't have any luck using the "::Rf_foo" syntax to access an R function from Rcpp code.

Here is a simplified version of my code below (yes, I am writing a gibbs sampler).

#include <Rcpp.h>                                                                                                                                     
using namespace Rcpp;                                                                                                                                 

// C++ implementation of the R which() function.                                                                                                      
int whichC(NumericVector x, double val) {                                                                                                             
  int ind = -1;                                                                                                                                       
  int n = x.size();                                                                                                                                   
  for (int i = 0; i < n; ++i) {                                                                                                                       
    if (x[i] == val) {                                                                                                                                
      if (ind == -1) {                                                                                                                                
        ind = i;                                                                                                                                      
      } else {                                                                                                                                        
        throw std::invalid_argument( "value appears multiple times." );                                                                               
      }                                                                                                                                               
    } // end if                                                                                                                                       
  } // end for                                                                                                                                        
  if (ind != -1) {                                                                                                                                    
    return ind;                                                                                                                                       
  } else {                                                                                                                                            
    throw std::invalid_argument( "value doesn't appear here!" );                                                                                      
    return -1;                                                                                                                                        
  }                                                                                                                                                   
}                                                                                                                                                     

// [[Rcpp::export]]                                                                                                                                   
int multSample(double p1, double p2, double p3) {                                                                                                     
  NumericVector params(3);                                                                                                                            
  params(0) = p1;                                                                                                                                     
  params(1) = p2;                                                                                                                                     
  params(2) = p3;                                                                                                                                     

  // HERE'S THE PROBLEM.                                                                                                                              
  RObject sampled = rmultinom(1, 1, params);                                                                                                          
  int out = whichC(as<NumericVector>(sampled), 1);                                                                                                    
  return out;                                                                                                                                         
}

I am new to c++, so I realize that alot of this code is probably nooby and inefficient. I'm open to suggestions on how to improve my c++ code, but my priority is hearing about the rmultinom business. Thanks!

Btw, I apologize for similarity with this thread, but

1. The answer didn't work for my purposes
2. The difference may be enough to warrant a different question (do you think so? )
3. The question was posted and answered a year ago.

Answer 1

The example and links provided by Kevin enabled me to find an answer that worked. There was some wrangling of types. I wrote a function that allows you to sample one vector from a multinomial distribution. The code is below.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector oneMultinomC(NumericVector probs) {
    int k = probs.size();
    SEXP ans;
    PROTECT(ans = RF_allocVector(INTSXP, k));
    probs = RF_coerceVector(probs, REALSXP);
    rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
    UNPROTECT(1);
    return ans;
}

I don't understand half of what is going on here. Particularly, I don't understand the fourth argument to 'rmultinom'. I know that it is a pointer to a memory location of where to store the output, but I don't understand the '[0]' bit. Nonetheless, it works. Gibbs sample away, boys and girls.

Answer 2

If I attempt to compile your code, I get the compiler error:

> Rcpp::sourceCpp('~/scratch/multSample.cpp')
multSample.cpp:33:21: error: no matching function for call to 'rmultinom'
  RObject sampled = rmultinom(1, 1, params);
                    ^~~~~~~~~
/Library/Frameworks/R.framework/Resources/include/Rmath.h:449:6: note: candidate function not viable: requires 4 arguments, but 3 were provided
void    rmultinom(int, double*, int, int*);
        ^
1 error generated.

As it suggests, you haven't specified the arguments correctly. Note that the rmultinom interface is a bit awkward compared to other functions: it fills the memory pointed at by *rn, rather than returning a new object (with its own freshly allocated memory).

If you look in the R sources you'll see the interface, and you can see an example of it being used here as well (in fact, stats makes a wrapper function that does some more argument checking and whatnot). But notice how it's used here:

rmultinom(size, REAL(prob), k, &INTEGER(ans)[ik]);

In other words, it's filling an INTSXP called ans by passing that address of that memory on to the rmultinom function.

So if you want to use this function from Rcpp, you'll have to do something similar -- but perhaps this deserves a similar sugar vectorization treatment to avoid the ugliness of that interface.

You might try doing something like:

IntegerMatrix sampled(nrow, ncol);
rmultinom(1, 1, params, sampled.begin());

or something to that effect.

Answer 3

Below are the answer by user95215 amended so that it compiles, and another version more in the Rcpp style:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector oneMultinomC(NumericVector probs) {
    int k = probs.size();
    SEXP ans;
    PROTECT(ans = Rf_allocVector(INTSXP, k));
    probs = Rf_coerceVector(probs, REALSXP);
    rmultinom(1, REAL(probs), k, &INTEGER(ans)[0]);
    UNPROTECT(1);
    return(ans);
}

// [[Rcpp::export]]
IntegerVector oneMultinomCalt(NumericVector probs) {
    int k = probs.size();
    IntegerVector ans(k);
    rmultinom(1, probs.begin(), k, ans.begin());
    return(ans);
}