RcppArmadillo's ability to fill matrix columns using rnorm()

Question

I'm trying to implement an rmvrnorm function that mimics the output from mvtnorm's rmvnorm function.

Here is what I have so far:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;

// [[Rcpp::export]]
arma::mat rmvrnorm_arma(int n, arma::vec mu, arma::mat sigma) {
   int ncols = sigma.n_cols;
   // works with rng?
   arma::mat Y(n,ncols);
   RNGScope scope;
   for(int i = 0; i<ncols; i++){
     // Fill a column using rnorm(n, mean = 0, sd = 1)
     Y.col(i) = rnorm(n);
   }

   return arma::repmat(mu, 1, n).t() + Y * arma::chol(sigma);
}

I am receiving an error on the line within the loop.

Specifically, the line:

         Y.col(i) = rnorm(n);

It states:

" rmvrnorm_arma.cpp:16:15: error: no viable overloaded '='

 Y.col(i) = rnorm(n);
 ~~~~~~~~ ^ ~~~~~~~~ "

I know arma has arma::randn() and I will potentially take a performance hit, however, using randn() negates R's RNG scope. I would like to be able to compare the output from this function with rmvnorm().

Answer 1

You need to help the compiler: use as to convert the default Rcpp::NumericVector object output from rnorm to an arma::vec object.

Y.col(i) = as<arma::vec>(rnorm(n));

Answer 2

Given how rnorm(int) works, i.e. it creates a NumericVector using the NormGenerator__mean0__sd1 class, which is just a wrapper around the ::norm_rand function, you could use something more STL-like:

std::generate( Y.col(i).begin(), Y.col(i).end(), ::norm_rand ) ;

Or you can directly use std::generate on the entire matrix.

std::generate( Y.begin(), Y.end(), ::norm_rand ) ;

This skips the creation (and memory allocation) of temporary NumericVector and arma::vec. You can easily derive from this approach if you need other parameters in rnorm by looking how various generators are implemented in Rcpp.

You could also use the imbue method in armadillo:

Y.imbue( norm_rand ) ;