Timing normal RNGs

January 16, 2013
By

(This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers)

In previous articles, we have seen that Rcpp can be particularly useful for simulations as it executes code at C++ speed. A very useful feature the API provided by R is the access to the R RNGs so that simulations at the C++ level can get precisely the same stream of random numbers as an R application would.

But sometimes that is not a requirement, and here will look into drawing normals from R, from the random number generator in Boost and the new one in C++11.

A first approach is by far the easiest: using Rcpp and its sugar function which reduces this to a single statement.

``````#include

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector rcppNormals(int n) {
RNGScope scope;             // also done by sourceCpp()
return rnorm(n);
}
``````

A quick test:

``````set.seed(42)
rcppNormals(10)
``````
``` [1]  1.37096 -0.56470  0.36313  0.63286  0.40427 -0.10612  1.51152
[8] -0.09466  2.01842 -0.06271
```

Next, the same via Boost. The caveats from the previous two Boost pieces apply: on some systems you may have to ensure access to the Boost headers, on some (such as my Linux system) it just works.

``````#include

#include
#include
#include

using namespace std;
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector boostNormals(int n) {

typedef boost::mt19937 RNGType; 	// select a generator, MT good default
RNGType rng(123456);			// instantiate and seed

boost::normal_distribution<> n01(0.0, 1.0);
boost::variate_generator< RNGType, boost::normal_distribution<> > rngNormal(rng, n01);

NumericVector V(n);
for ( int i = 0; i < n; i++ ) {
V[i] = rngNormal();
};

return V;
}
``````

A second test:

``````boostNormals(10)
``````
``` [1]  0.8400  0.8610  2.0907 -0.4437 -0.1029  1.5609  1.3874 -1.0453
[9] -1.6558  1.6198
```

And now for the same using the random number generator added to C++11. Here, the same caveats apply as before: we need to enable the C++11 extensions:

``````Sys.setenv("PKG_CXXFLAGS"="-std=c++11")
# or
# Sys.setenv("PKG_CXXFLAGS"="-std=c++0x")
``````

That way, we can compile this code:

``````#include
#include

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector cxx11Normals(int n) {

std::mt19937 engine(42);
std::normal_distribution<> normal(0.0, 1.0);

NumericVector V(n);
for ( int i = 0; i < n; i++ ) {
V[i] = normal(engine);
};

return V;
}
``````

And run it:

``````cxx11Normals(10)
``````
``` [1] -0.5502  0.5154  0.4739  1.3685 -0.9168 -0.1241 -2.0110 -0.4928
[9]  0.3926 -0.9292
```

Lastly, we can compare the runtime for these three in a quick benchmark study.

``````library(rbenchmark)

n <- 1e5

res <- benchmark(rcppNormals(n),
boostNormals(n),
cxx11Normals(n),
order="relative",
replications = 500)
print(res[,1:4])
``````
```             test replications elapsed relative
3 cxx11Normals(n)          500   3.764    1.000
1  rcppNormals(n)          500   4.102    1.090
2 boostNormals(n)          500   4.353    1.156
```

In this particularly example, all the RNGs take roughly the same time. It would be interesting to see how the Ziggurat algorithm (which is known to produce Normals rather fast) would fare.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...