The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Flexible and consistent simulation of a matrix of Monte Carlo variates

Riccardo Porreca, Roland Schmid

2022-03-14

Consider the Monte Carlo simulation of a matrix of i.i.d. normal random variables. We will show how rTRNG can be used to perform a consistent (fair-playing) simulation of a subset of the variables and simulations.

Consistent simulation in R

We rely on the TRNG engines exposed to R as reference classes by rTRNG.

library(rTRNG)

The mcMatR function below performs the full sequential Monte Carlo simulation of nrow normal i.i.d. samples of ncol variables using the yarn2 generator.

mcMatR <- function(nrow, ncol) {
  r <- yarn2$new(12358)
  M <- matrix(rnorm_trng(nrow * ncol, engine = r),
              nrow = nrow, ncol = ncol, byrow = TRUE)
  M
}

A second function mcSubMatR relies on jump and split operations to perform only a chunk [startRow, endRow] of simulations for a subset subCols of the variables.

mcSubMatR <- function(nrow, ncol,
                      startRow, endRow, subCols) {
  r <- yarn2$new(12358)
  r$jump((startRow - 1)*ncol)
  nSubCols <- endRow - startRow + 1
  S <- matrix(0.0, nrow, ncol)
  S[startRow:endRow, subCols] <-
    vapply(subCols,
           function(j) {
             rj = r$copy()
             rj$split(ncol, j)
             rnorm_trng(nSubCols, engine = rj)
           },
           FUN.VALUE = numeric(nSubCols))
  S
}

The parallel nature of the yarn2 generator ensures the sub-simulation obtained via mcSubMatR is consistent with the full sequential simulation.

rows <- 9
cols <- 5
M <- mcMatR(rows, cols)
startRow <- 4
endRow <- 6
subCols <- c(2, 4:5)
S <- mcSubMatR(rows, cols,
               startRow, endRow, subCols)
identical(M[startRow:endRow, subCols],
          S[startRow:endRow, subCols])
## [1] TRUE
M.1 M.2 M.3 M.4 M.5 S.1 S.2 S.3 S.4 S.5
1 0.20256 -0.41401 -0.76749 -0.33344 0.10718 0 0.00000 0 0.00000 0.00000
2 -2.76439 -1.15524 -0.39394 -1.16604 1.61759 0 0.00000 0 0.00000 0.00000
3 -0.42199 -1.13148 -0.30448 0.12741 -0.16111 0 0.00000 0 0.00000 0.00000
4 -0.94448 -1.86384 -1.03244 -0.41155 1.31009 0 -1.86384 0 -0.41155 1.31009
5 -0.09614 -0.16366 -0.31964 0.87053 0.77996 0 -0.16366 0 0.87053 0.77996
6 1.42049 -0.73062 -1.19459 -1.02146 0.07202 0 -0.73062 0 -1.02146 0.07202
7 -0.61202 0.02906 -0.29100 0.10095 -0.74647 0 0.00000 0 0.00000 0.00000
8 1.10246 -0.50507 0.01286 0.63140 -1.28893 0 0.00000 0 0.00000 0.00000
9 -0.08732 -0.32545 0.29099 0.62003 -0.94617 0 0.00000 0 0.00000 0.00000

Consistent simulation with Rcpp

We now use Rcpp to define functions mcMatRcpp and mcSubMatRcpp for the full sequential simulation and the sub-simulation, respectively. The Rcpp::depends attribute makes sure the TRNG library and headers shipped with rTRNG are available to the C++ code. Moreover, Rcpp::plugins(cpp11) enforces the C++11 standard required by TRNG >= 4.22.

// [[Rcpp::depends(rTRNG)]]
// TRNG >= 4.22 requires C++11
// [[Rcpp::plugins(cpp11)]]
#include <Rcpp.h>
#include <trng/normal_dist.hpp>
#include <trng/yarn2.hpp>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix mcMatRcpp(const int nrow, const int ncol) {
  NumericMatrix M(nrow, ncol);
  trng::yarn2 r(12358);
  trng::normal_dist<> normal(0.0, 1.0);
  for (int i = 0; i < nrow; i++) {
    for (int j = 0; j < ncol; j++) {
      M(i, j) = normal(r);
    }
  }
  return M;
}
// [[Rcpp::export]]
NumericMatrix mcSubMatRcpp(const int nrow, const int ncol,
                           const int startRow,
                           const int endRow,
                           const IntegerVector subCols) {
  NumericMatrix M(nrow, ncol);
  trng::yarn2 r(12358), rj;
  trng::normal_dist<> normal(0.0, 1.0);
  r.jump((startRow - 1) * ncol);
  for (IntegerVector::const_iterator jSub = subCols.begin();
       jSub < subCols.end(); jSub++) {
    int j = *jSub - 1;
    rj = r;
    rj.split(ncol, j);
    for (int i = startRow - 1; i < endRow; i++) {
      M(i, j) = normal(rj);
    }
  }
  return M;
}

As seen above for the R case, consistency of the simulation obtained via mcSubMatRcpp with the full sequential simulation is guaranteed.

rows <- 9
cols <- 5
startRow <- 4
endRow <- 6
subCols <- c(2, 4:5)
M <- mcMatRcpp(rows, cols)
S <- mcSubMatRcpp(rows, cols, startRow, endRow, subCols)
identical(M[startRow:endRow, subCols],
          S[startRow:endRow, subCols])
## [1] TRUE
M.1 M.2 M.3 M.4 M.5 S.1 S.2 S.3 S.4 S.5
1 0.20256 -0.41401 -0.76749 -0.33344 0.10718 0 0.00000 0 0.00000 0.00000
2 -2.76439 -1.15524 -0.39394 -1.16604 1.61759 0 0.00000 0 0.00000 0.00000
3 -0.42199 -1.13148 -0.30448 0.12741 -0.16111 0 0.00000 0 0.00000 0.00000
4 -0.94448 -1.86384 -1.03244 -0.41155 1.31009 0 -1.86384 0 -0.41155 1.31009
5 -0.09614 -0.16366 -0.31964 0.87053 0.77996 0 -0.16366 0 0.87053 0.77996
6 1.42049 -0.73062 -1.19459 -1.02146 0.07202 0 -0.73062 0 -1.02146 0.07202
7 -0.61202 0.02906 -0.29100 0.10095 -0.74647 0 0.00000 0 0.00000 0.00000
8 1.10246 -0.50507 0.01286 0.63140 -1.28893 0 0.00000 0 0.00000 0.00000
9 -0.08732 -0.32545 0.29099 0.62003 -0.94617 0 0.00000 0 0.00000 0.00000

Consistent parallel simulation with RcppParallel

The same technique used for generating a sub-set of the simulations can be exploited for performing a parallel simulation in C++. We can embed the body of mcSubMatRcpp above into an RcppParallel::Worker for performing chunks of Monte Carlo simulations in parallel, for any subset subCols of the variables.

struct MCMatWorker : public Worker
{
  RMatrix<double> M;
  const RVector<int> subCols;

  // constructor
  MCMatWorker(NumericMatrix M,
              const IntegerVector subCols)
    : M(M), subCols(subCols) {}

  // operator processing an exclusive range of indices
  void operator()(std::size_t begin, std::size_t end) {
    trng::yarn2 r(12358), rj;
    trng::normal_dist<> normal(0.0, 1.0);
    r.jump((int)begin * M.ncol());
    for (IntegerVector::const_iterator jSub = subCols.begin();
         jSub < subCols.end(); jSub++) {
      int j = *jSub - 1;
      rj = r;
      rj.split(M.ncol(), j);
      for (int i = (int)begin; i < (int)end; i++) {
        M(i, j) = normal(rj);
      }
    }
  }
};
// [[Rcpp::export]]
NumericMatrix mcMatRcppParallel(const int nrow, const int ncol,
                                const IntegerVector subCols) {
  NumericMatrix M(nrow, ncol);
  MCMatWorker w(M, subCols);
  parallelFor(0, M.nrow(), w);
  return M;
}

The parallel nature of the yarn2 generator ensures the parallel simulation is playing fair, i.e. is consistent with the sequential simulation.

M <- mcMatRcpp(rows, cols)
Mp <- mcMatRcppParallel(rows, cols, seq_len(ncol(M)))
identical(M, Mp)
## [1] TRUE

Similarly, we can achieve a consistent parallel simulation of a subset of the variables only.

Sp <- mcMatRcppParallel(rows, cols, subCols)
identical(M[, subCols],
          Sp[, subCols])
## [1] TRUE
M.1 M.2 M.3 M.4 M.5 Sp.1 Sp.2 Sp.3 Sp.4 Sp.5
1 0.20256 -0.41401 -0.76749 -0.33344 0.10718 0 -0.41401 0 -0.33344 0.10718
2 -2.76439 -1.15524 -0.39394 -1.16604 1.61759 0 -1.15524 0 -1.16604 1.61759
3 -0.42199 -1.13148 -0.30448 0.12741 -0.16111 0 -1.13148 0 0.12741 -0.16111
4 -0.94448 -1.86384 -1.03244 -0.41155 1.31009 0 -1.86384 0 -0.41155 1.31009
5 -0.09614 -0.16366 -0.31964 0.87053 0.77996 0 -0.16366 0 0.87053 0.77996
6 1.42049 -0.73062 -1.19459 -1.02146 0.07202 0 -0.73062 0 -1.02146 0.07202
7 -0.61202 0.02906 -0.29100 0.10095 -0.74647 0 0.02906 0 0.10095 -0.74647
8 1.10246 -0.50507 0.01286 0.63140 -1.28893 0 -0.50507 0 0.63140 -1.28893
9 -0.08732 -0.32545 0.29099 0.62003 -0.94617 0 -0.32545 0 0.62003 -0.94617

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.