Search code examples
c++rrcppbayesianmcmc

Using rmultinom with Rcpp


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.

Solution

  • 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);
    }