It is well-known that if you are writing simulation code in R you can often gain a performance boost by rewriting parts of your simulation in C++. These days the easiest way to do that of course is to use Rcpp. Simulation usually depends on random variates, and usually great numbers of them. One of the issues that may arise is that your simulation needs to execute on the C++ side of things. For example, if you decide to programme your Metropolis-Hastings algorithm (not technically a simulation I know) in Rcpp, then you are going to need to be able to generate hundreds of thousands, if not millions, of random numbers. You can use Rcpp’s features to call R routines from within Rcpp to do this, e.g.
Function rnorm("rnorm"); rnorm(100, _["mean"] = 10.2, _["sd"] = 3.2 );
but this has a certain overhead. C++ has had built-in in random number generation functionality since at least the C+11 standard (and probably since the C+0X standard). The random
header file provides a Mersenne-Twister uniform random number generator (RNG), a Linear Congruential Generator (LCG), and a Subtract-with-Carry RNG. There is also a variety of standard distributions available, described here.
Uniform random variates
The ability to generate good quality uniform random variates is essential, and the mt19937
engine provides. The 19937 refers to the Mersenne Prime \((2^{19937}-1)\) that this algorithm is based on, and also to its period length. There are four steps required to generate uniform random variates. These are:
- Include the
random
header file - Construct an
mt19937
random number engine, and initialise it with a seed - Construct a \(U(0,1)\) random number generator
- Use your engine and your uniform random number generator to draw variates
In code we would write
#include <random> #include <Rcpp.h> using namespace std; using namespace Rcpp; mt19937 mtEngine; uniform_real_distribution<double> rngU; //[[Rcpp::export]] void setSeed(unsigned int seed){ mtEngine = mt19937(seed); rngU = uniform_real_distribution<>(0.0, 1.0); } double runif(void){ return rngU(mtEngine); }
The function runif
can now be called with runif(). Note that the setSet
function has been exported so that you can initialize the RNG engine with a seed of your choice.
How about normal random variates?
It does not require very much more effort to add a normal RNG to your code. We simply add
normal_distribution<double> rngZ;
to our declared variables, and
//[[Rcpp::export]] void setSeed(unsigned int seed){ mtEngine = mt19937(seed); rngU = uniform_real_distribution<>(0.0, 1.0); rngZ = normal_distribution<double>(0.0, 1.0); } double rnorm(double mu = 0, double sigma = 1){ return rngZ(mtEngine) * sigma + mu; }
to our code base. Now rnorm
can be called without arguments to get standard (\(N(0,1)\)) random variates, or with a mean, or a standard deviation, or both to get \(N(\mu,\sigma^2)\) random variates
Rcpp does it
No doubt someone is going to tell me that Romain and Dirk have thought of this already for you, and that my solution is unnecessary Morris Dancing. However, I think there is merit in knowing how to use the standard C++ libraries.
Please note that I do not usually advocate having global variables such as those in the code above. I would normally make mtEngine
, rngU
, and rngZ
private member variables a class and then either instantiate it using an exported Rcpp function, or export the class and essential functions using an Rcpp module.
Working C++ code and an R test script can be found here in the RNG
folder. Enjoy!