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.
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.
We rely on the TRNG engines exposed to R as reference classes by 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 |
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.
// [[Rcpp::depends(rTRNG)]]
#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 |
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.
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.