通过RcppEigen随机状态设置随机

Question

有没有办法通过RcppEigen将随机状态传递给Eigen的setRandom？还是我需要使用runif？

下面是一个例子：

// [[Rcpp::depends(RcppEigen)]] 
#include <RcppEigen.h> 
using namespace Rcpp; 
using Eigen::MatrixXd; 
using Eigen::VectorXd; 

// [[Rcpp::export]] 
NumericVector fx() { 
    RNGScope scope; 
    MatrixXd x(3,2); 
    x=x.setRandom(); 
    x.col(1)=as<VectorXd>(runif(3,0,1)); 

return wrap(x); 
} 

测试它：

set.seed(42); fx() 
#   [,1]  [,2] 
#[1,] -0.8105760 0.9148060 
#[2,] 0.6498853 0.9370754 
#[3,] 0.6221027 0.2861395 

set.seed(42); fx() 
#   [,1]  [,2] 
#[1,] -0.9449154 0.9148060 
#[2,] 0.8063267 0.9370754 
#[3,] -0.0673205 0.2861395 

注如何塔2（即，runif）是可再现的，但第1列（即，setRandom）不是。

Answer 1

Eigen中的RNG与R的RNG正交。如你所知，我们有RNGScope来处理R，并允许你set.seed();你对Eigen的RNG做了什么是分开的。

特别是，您需要添加胶水以将Eigen的RNG注册为R的用户提供的RNG。默认情况下，R不知道Eigen，这对于R用户期望R的RNG是有意义的。

你的问题在这里似乎是一个vanilla'Eigen编程'问题，你无法初始化特征RNG。它与Rcpp无关。

编辑：这是你的例子，更正。原来Eigen使用系统RNG，因此您需要srand()来播种它。响应

// [[Rcpp::depends(RcppEigen)]] 
#include <RcppEigen.h> 
using namespace Rcpp; 
using Eigen::MatrixXd; 
using Eigen::VectorXd; 

// [[Rcpp::export]] 

    RNGScope scope; 
    MatrixXd x(3,2); 
    srand(seed); 
    x=x.setRandom(); 
    x.col(1)=as<VectorXd>(runif(3,0,1)); 

    return wrap(x); 
} 

编辑2 OP的评论：

R> sourceCpp("/tmp/roland.cpp") 
R> set.seed(42); fx(42) 
      [,1]  [,2] 
[1,] -0.933060 0.914806 
[2,] -0.340072 0.937075 
[3,] 0.381271 0.286140 
R> set.seed(42); fx(42) 
      [,1]  [,2] 
[1,] -0.933060 0.914806 
[2,] -0.340072 0.937075 
[3,] 0.381271 0.286140 
R> 

它使用你的代码这款改装版

我不希望Rcpp::runif()和太大的速度差srand() Eigen所做的调用，所以你仍然坚持srand()一直有问题，并可能在系统之间有不同的表现。

快速演示脚本：

#include <RcppEigen.h> 

// [[Rcpp::depends(RcppEigen)]] 

// [[Rcpp::export]] 
Rcpp::NumericVector v1(int n) { 
    return Rcpp::runif(n); 
} 

// [[Rcpp::export]] 
Rcpp::NumericVector v2(int n) { 
    Eigen::VectorXd x(n); 
    x = x.setRandom(); 
    return Rcpp::wrap(x); 
} 

/*** R 
library(rbenchmark) 
N <- 1e7 
benchmark(v1(N), v2(N)) 
*/ 

产生

R> sourceCpp("/tmp/roland.cpp") 

R> library(rbenchmark) 

R> N <- 1e7 

R> benchmark(v1(N), v2(N)) 
    test replications elapsed relative user.self sys.self user.child sys.child 
1 v1(N)   100 12.633 1.000 11.356 1.261   0   0 
2 v2(N)   100 17.222 1.363 13.981 3.198   0   0 
R> 

，并注意RcppEigen是慢这里即使在简单设置只是创造载体。但我们在这里说的是微秒，这可能不是我会担心的任何方式在真正的应用程序，这很可能有其他瓶颈。