Passing a vector of lambdas to Rcpp's rpois

Question

Within R, rpois can be passed a vector of lambdas describing multiple Poisson distributions, e.g.

rpois(5, (1:5)*1000)

# [1] 1043 1974 3002 3930 4992

Above, each element of the output vector is drawn from a different Poisson distribution, with means of 1000, 2000, 3000, 4000 and 5000, respectively.

If I have an arma::mat (using these because elsewhere I am using cubes) containing lambdas of Poisson distributions, what is the best way to pass these (one row at a time) to rpois within Rcpp?

Here's a toy example, and an excerpt from the ensuing error message:

library(inline)
library(RcppArmadillo)
code <- "
  using namespace Rcpp;
  using namespace arma;
  arma_rng::set_seed(42); // Dirk's seed of choice

  mat lam = randu(5, 5); // ignore the fact these are all 0-1
  mat out(5, 5);

  for (int i = 0; i < 5; i++) {
    out.row(i) = rpois(5, lam.row(i));
  }      

  return(wrap(out));
"

f <- cxxfunction(body=code, plugin="RcppArmadillo")

# cannot convert 'arma::subview_row<double>' to 'double' for argument '2' 
#   to 'Rcpp::NumericVector Rcpp::rpois(int, double)'

I must admit my understanding of type conversion in c++ is quite poor. Is what I'm trying to do possible (my guess is no, since it seems rpois is expecting a double), or do I need to iterate over each cell of the matrix, generating a single deviate each time?

Answer 1

From C/C++, you have access to at least 2 Poisson deviate generation routines (search for rpois in this online manual).

Their declarations are as follows:

double R::rpois(double mu);
NumericVector Rcpp::rpois(int n, double mu);

None of them allows for passing > 1 values in mu (a.k.a. lambda). The first function is the R's original routine, used to implement rpois as we know from the R stats package (the one that is vectorized w.r.t. all its arguments). Given a single mu, it returns a single (pseudo)random deviate.

The second one is the so-called Rcpp sugar function. It allows to calculate n deviates at a time and return them as a NumericVector (again, by using R::rpois).

In other words, you should rather use two nested for loops to fill a matrix, by calling R::rpois. Don't be afraid of such an approach, this is C++. :)