将 dnorm 和 pnorm 函数从 R 导入到 Rcpp
Issue importing dnorm and pnorm functions from R to Rcpp
我在尝试将 dnorm()
和 pnorm()
函数从 R 包 stats
导入 Rcpp 时遇到了以下奇怪的行为。
我将均值 0 和标准差 1 的 dnorm()
和 pnorm()
应用于观察向量。如果我对同一个向量重复多次计算,有时会得到不同的结果。看起来我在导入这两个函数时做错了,但这并不能解释为什么结果有时好,而其他不一致。
我把从 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
}
有时,实验在到达循环末尾之前就中断了。有时它没有。没有明确的模式。
感谢您的回复。
天哪,这些已经可用了。
首先,通过糖:
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
,当然参数不同。)
我在尝试将 dnorm()
和 pnorm()
函数从 R 包 stats
导入 Rcpp 时遇到了以下奇怪的行为。
我将均值 0 和标准差 1 的 dnorm()
和 pnorm()
应用于观察向量。如果我对同一个向量重复多次计算,有时会得到不同的结果。看起来我在导入这两个函数时做错了,但这并不能解释为什么结果有时好,而其他不一致。
我把从 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
}
有时,实验在到达循环末尾之前就中断了。有时它没有。没有明确的模式。
感谢您的回复。
天哪,这些已经可用了。
首先,通过糖:
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
,当然参数不同。)