dqrng с 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 .Если я правильно понял документацию, я могу заменить R::rnorm(2.0, 1.0) на dqrng::dqrnorm(1, 2.0, 1.0)[1].Это правильно?Как насчет замены R::rbinom(1,out(row, 3))?Мне не удалось найти ссылки на рисунок из бинома в документации

---

Мне удалось написать следующий код, который рисует из биномиального распределения параллельно:

#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 объект и перенастроить его в цикле, см. здесь для аналогичного вопроса.