c++ - 将 dnorm 和 pnorm 函数从 R 导入到 Rcpp

Question

在尝试将 R 包 stats 中的 dnorm() 和 pnorm() 函数导入 Rcpp 时，我遇到了以下奇怪的行为。我将均值为 0 且标准差为 1 的 dnorm() 和 pnorm() 应用于观察 vector 。如果我对同一个 vector 重复多次计算，有时会得到不同的结果。看起来我在两个函数的导入上做错了，但这并不能解释为什么结果有时很好，而其他不一致。

我把从 stats 导入 dnorm() 和 pnorm() 函数的函数放在这里:

    {test_imp_script<-'
     using Eigen::VectorXd;
     VectorXd a = Rcpp::as<VectorXd>(aa);

     Environment st("package:stats");
     Function dn = st["dnorm"];
     Function pn = st["pnorm"]; 
     SEXP dd_a = dn(a,0,1);
     SEXP pp_a = pn(a,0,1);
     VectorXd d_a = as<VectorXd>(dd_a);
     VectorXd p_a = as<VectorXd>(pp_a);

     return List::create(Named("d") = d_a,Named("p") = p_a);'}
    test_imp <- cxxfunction(signature( aa="numeric"), test_imp_script, plugin = "RcppEigen")

现在，尝试使用 2000 次 test_imp 重复运行以下示例。

    set.seed(123)
    t<-rnorm(1000,0,1)
    a<-test_imp(t)
    for (i in 1:2000){
      set.seed(123)
      b<-test_imp(t)
      cat(i,"d",c(sum(b$d),sum(b$p)),"\n")
      if (any(a$d!=b$d)) break    
      if (any(a$p!=b$p)) break    
    }

有时，实验在到达循环末尾之前就中断了。有时它没有。没有明确的模式。

感谢您的回复。

Answer 1

天哪，这些已经可用了。

首先，通过糖:

R> library(Rcpp)
R> cppFunction("NumericVector ex1() { return rnorm(5, 0, 1); }")
R> set.seed(42); ex1()
[1]  1.370958 -0.564698  0.363128  0.632863  0.404268
R> 

其次，通过方便的命名空间(但请注意，这不是矢量化的)

R> cppFunction("double ex2() { return R::rnorm(0, 1); }")
R> set.seed(42); ex2()
[1] 1.37096
R> 

第三，如果你坚持，直接调用R API 越少越难看之前的答案:

R> cppFunction("double ex3() { return Rf_rnorm(0, 1); }")
R> set.seed(42); ex3()
[1] 1.37096
R> 

我们竭尽全力提供免费且全面的文档。我真的不认为你看过它。这样做对您有利，因为您将不再依赖其他人在这里为您重新输入。

(是的，您确实要求了 dnorm 和 pnorm 但 all 所谓的 d/p/q/r 统计分布的函数是并行实现的。因此，在 rnorm 存在的地方，您也有 dnorm、qnorm 和 pnorm ，当然有不同的参数。)