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

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

我使用以下代码编写了R::rnorm and R::rbinom我理解这是一个不要那样做 https://stackoverflow.com/a/54717911/1800784.

#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 https://cran.r-project.org/web/packages/dqrng/vignettes/dqrng.html。如果我正确理解文档，我可以替换R::rnorm(2.0, 1.0) with dqrng::dqrnorm(1, 2.0, 1.0)[1]。是对的吗？更换怎么样R::rbinom(1,out(row, 3))？我无法找到任何关于二项式绘图的参考文档 https://cran.r-project.org/web/packages/dqrng/dqrng.pdf

我能够编写以下并行从二项式分布中提取的代码：

#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 的函数？

我写的代码有什么问题吗？

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

// [[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 and j。保留一个全局可能会更有效dist对象并在循环内重新配置它，请参阅here https://stackoverflow.com/questions/14268049/how-to-set-the-parameters-of-boostuniform-int-distribution-outside-of-the-stan对于类似的问题。