Rcpp - генерирует несколько случайных наблюдений из пользовательского распределения - PullRequest
Новая
       5

Rcpp - генерирует несколько случайных наблюдений из пользовательского распределения

0 голосов
/ 07 сентября 2018

Этот вопрос относится к предыдущему при вызове функций внутри функций в Rcpp.

Мне нужно сгенерировать большое количество случайных отрисовок из пользовательского дистрибутива способом, подобным rnorm () или rbinom (), с дополнительным усложнением, что моя функция выдает векторный вывод.

В качестве решения я подумал об определении функции, которая генерирует наблюдения из пользовательского распределения, а затем о главной функции, которая извлекает n раз из генерирующей функции через цикл for. Ниже приведена значительно упрощенная рабочая версия кода: 

#include <Rcpp.h>
using namespace Rcpp;

// generating function
NumericVector gen(NumericVector A, NumericVector B){
  NumericVector out = no_init_vector(2); 
  out[0] = R::runif(A[0],A[1]) + R::runif(B[0],B[1]);
  out[1] = R::runif(A[0],A[1]) - R::runif(B[0],B[1]);
  return out;
}

// [[Rcpp::export]]
// draw n observations
NumericVector rdraw(int n, NumericVector A, NumericVector B){
  NumericMatrix out = no_init_matrix(n, 2);
  for (int i = 0; i < n; ++i) {
    out(i,_) = gen(A, B); 
  }
  return out;
}

Я ищу способы ускорить розыгрыш. Мои вопросы: есть ли более эффективная альтернатива циклу for? Поможет ли распараллеливание в этом случае?

Спасибо за любую помощь!

1 Ответ

0 голосов
/ 07 сентября 2018

Есть несколько способов ускорить это:

  1. Используйте inline на gen(), уменьшая количество вызовов функций.
  2. Используйте Rcpp::runif вместо цикла с R::runif для удаления еще большего количества вызовов функций.
  3. Используйте более быстрый ГСЧ, который допускает параллельное выполнение.

Здесь пункты 1. и 2.:

#include <Rcpp.h>
using namespace Rcpp;

// generating function
inline NumericVector gen(NumericVector A, NumericVector B){
  NumericVector out = no_init_vector(2); 
  out[0] = R::runif(A[0],A[1]) + R::runif(B[0],B[1]);
  out[1] = R::runif(A[0],A[1]) - R::runif(B[0],B[1]);
  return out;
}

// [[Rcpp::export]]
// draw n observations
NumericVector rdraw(int n, NumericVector A, NumericVector B){
  NumericMatrix out = no_init_matrix(n, 2);
  for (int i = 0; i < n; ++i) {
    out(i,_) = gen(A, B); 
  }
  return out;
}

// [[Rcpp::export]]
// draw n observations
NumericVector rdraw2(int n, NumericVector A, NumericVector B){
  NumericMatrix out = no_init_matrix(n, 2);
  out(_, 0) = Rcpp::runif(n, A[0],A[1]) + Rcpp::runif(n, B[0],B[1]);
  out(_, 1) = Rcpp::runif(n, A[0],A[1]) - Rcpp::runif(n, B[0],B[1]);
  return out;
}

/*** R
set.seed(42)
system.time(rdraw(1e7, c(0,2), c(1,3)))
system.time(rdraw2(1e7, c(0,2), c(1,3)))
*/

Результат: 

> set.seed(42)

> system.time(rdraw(1e7, c(0,2), c(1,3)))
   user  system elapsed 
  1.576   0.034   1.610 

> system.time(rdraw2(1e7, c(0,2), c(1,3)))
   user  system elapsed 
  0.458   0.139   0.598

Для сравнения, ваш исходный код занял около 1,8 с для 10 ^ 7 розыгрышей. Для пункта 3. Я адаптирую код из параллельной виньетки моего dqrng пакета: 

#include <Rcpp.h>
// [[Rcpp::depends(dqrng)]]
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
Rcpp::NumericMatrix rdraw3(int n, Rcpp::NumericVector A, Rcpp::NumericVector B, int seed, int ncores) {
  dqrng::uniform_distribution distA(A(0), A(1));
  dqrng::uniform_distribution distB(B(0), B(1));
  dqrng::xoshiro256plus rng(seed);
  Rcpp::NumericMatrix res = Rcpp::no_init_matrix(n, 2);
  RcppParallel::RMatrix<double> output(res);

  #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 genA = std::bind(distA, std::ref(lrng));
  auto genB = std::bind(distB, std::ref(lrng));      

  #pragma omp for
  for (int i = 0; i < n; ++i) {
    output(i, 0) = genA() + genB();
    output(i, 1) = genA() - genB();
  }
  }
  return res;
}

/*** R
system.time(rdraw3(1e7, c(0,2), c(1,3), 42, 2))
*/

Результат: 

> system.time(rdraw3(1e7, c(0,2), c(1,3), 42, 2))
   user  system elapsed 
  0.276   0.025   0.151

Таким образом, с более быстрым ГСЧ и умеренным параллелизмом мы можем получить порядок времени выполнения. Конечно, результаты будут другими, но итоговая статистика должна быть одинаковой.

...