问题描述
我正在尝试学习如何从 Rcpp OpenMP 循环中的正态分布和二项分布中抽取随机数.
I'm trying to learn how to draw random numbers from a normal and a binomial distribution from within a Rcpp OpenMP loop.
我使用 R::rnorm 和 R::rbinom 编写了以下代码，我认为这是一个 不要那样做.
I wrote the following code using R::rnorm and R::rbinom which I understand to be a be a don't do that.
#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)) 怎么样?我在 文档中找不到任何关于从二项式绘图的参考
As suggested by @ralf-stubner I'm trying to replace that code with 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
我能够编写以下从二项式分布中并行抽取的代码:
I was able to write the following code that draws from a binomial distribution in parallel:
#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的函数??
Is there a way to call out(j,i) = gen(); and make the probability be a function of i and j??
我写的代码有问题吗?
推荐答案
一个简单的解决方案是在循环内创建一个新的分发对象:
A simple solution would be to create a new distribution object within the loop:
// [[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 对象并在循环中重新配置它可能更有效，请参阅 此处 类似问题.
This way you can make p depend on i and j. It might be more efficient to keep one global dist object and reconfigure it within the loop, see here for a similar question.
这篇关于dqrng 与 Rcpp 用于绘制正态分布和二项分布的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！