dqrng with Rcpp 用于从正态分布和二项分布中绘制

Question

我正在尝试学习如何从 Rcpp OpenMP 循环中的正态分布和二项分布中提取随机数。

我使用 R::rnorm 和 R::rbinom 编写了以下代码，据我所知是 .

#include <RcppArmadillo.h>
#include <omp.h>
#include <dqrng.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(dqrng)]]
// [[Rcpp::plugins(openmp)]]
using namespace Rcpp;

// [[Rcpp::export]]
arma::mat  my_matrix3(const arma::mat& A, 
                      const arma::mat& B) {

  dqrng::dqRNGkind("Xoroshiro128+");
  dqrng::dqset_seed(42);



  const int nObservations = A.n_cols;
  const int nDraws = B.n_rows;
  const int nRows = nObservations * nDraws;

  // Show initialization information
  Rcpp::Rcout << "nObservations: " << nObservations << std::endl 
              << "nDraws: " << nDraws << std::endl 
              << "nRows: " << nRows << std::endl;

  arma::mat out(nRows, 5);

  // Show trace of matrix construction
  Rcpp::Rcout << "out - rows: " << out.n_rows << std::endl 
              << "out - columns: " << out.n_cols << std::endl;

  int i,n,iter,row;
  double dot_product, random_number, p;
  omp_set_num_threads(2);
#pragma omp parallel for private(i, n, iter, row)
  for(i = 0; i < nDraws; ++i){
    for(n = 0; n < nObservations; ++n) {
      row = i * nObservations + n;
      dot_product = arma::as_scalar(A.col(n).t() * B.row(i).t());
      // Show trace statement of index being accessed
      out(row, 0) = i + 1;
      out(row, 1) = n + 1;
      out(row, 2) = R::rnorm(2.0, 1.0) ;//dqrng::dqrnorm(1, 2.0, 1.0)[1]; 
      out(row, 3) = 1 / (1 + std::exp(-out(row, 3) - std::exp(dot_product)));
      out(row, 4) = R::rbinom(1,out(row, 3));
    }
  }

  return out;
}

/*** R
set.seed(9782)
A <- matrix(rnorm(10), nrow = 5)
B <- matrix(rnorm(10), ncol = 5)
test <- my_matrix3(A = A, B = B)
test

mean(test[,5])
*/

按照@ralf-stubner 的建议，我正在尝试用 dqrng. If I understood the documentation correctly, I can replace R::rnorm(2.0, 1.0) with dqrng::dqrnorm(1, 2.0, 1.0)[1]. Is that right? What about replacing R::rbinom(1,out(row, 3))? I was not able to find any reference to drawing from a binomial in the documentation

替换该代码

我能够编写以下并行二项分布的代码：

#include <RcppArmadillo.h>
// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <pcg_random.hpp>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
  //dqrng::uniform_distribution dist(0.0, 1.0);
  boost::random::binomial_distribution<int> dist(1,p);
  dqrng::xoshiro256plus rng(42);
  arma::mat out(n,m);
  // ok to use rng here



#pragma omp parallel num_threads(ncores)
{
  dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
  lrng.jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps 
  auto gen = std::bind(dist, std::ref(lrng));

#pragma omp for
  for (int i = 0; i < m; ++i) {
    double lres(0);
    for (int j = 0; j < n; ++j) {
      out(j,i) = gen(); /// CAN I MAKE P BE A FUNCTION OF i and j? how???
    }
  }
}
// ok to use rng here
return out;
}

/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/

> parallel_random_matrix(5, 5, 4, 0.75)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1
[2,]    0    1    1    1    1
[3,]    0    1    0    1    0
[4,]    1    1    1    1    1
[5,]    1    1    1    1    1

有没有办法调用 out(j,i) = gen(); 并使概率成为 i 和 j 的函数？

我写的代码有问题吗？

Answer 1

一个简单的解决方案是在循环内创建一个新的分布对象：

// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
  dqrng::xoshiro256plus rng(42);
  arma::mat out(n,m);
  // ok to use rng here

#pragma omp parallel num_threads(ncores)
{
  dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
  lrng.jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps 

#pragma omp for
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      // p can be a function of i and j
      boost::random::binomial_distribution<int> dist(1,p);
      auto gen = std::bind(dist, std::ref(lrng));
      out(j,i) = gen();
    }
  }
}
  // ok to use rng here
  return out;
}

/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/

这样你就可以让 p 依赖于 i 和 j。保留一个全局 dist 对象并在循环中重新配置它可能更有效，请参阅 here 了解类似问题。

dqrng with Rcpp 用于从正态分布和二项分布中绘制

dqrng with Rcpp for drawing from a normal and a binomial distribution

rcpp