通過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是慢這裏即使在簡單設置只是創造載體。但我們在這裏說的是微秒，這可能不是我會擔心的任何方式在真正的應用程序，這很可能有其他瓶頸。