c++ - 计算 vector 正态分布的cdf的最快方法-R::pnorm vs erfc vs？

我希望我改写的问题现在符合Stackoverflow的标准。请考虑以下示例。我正在编写一个Log-Likelihood函数，其中在 vector 上计算cdf是最耗时的部分。示例1使用R::pnorm，示例2使用erfc近似普通cdf。如您所见，结果足够相似，ercf版本要快一些。

然而，在实践中(在MLE内)事实证明ercf并不那么精确，这使得该算法可以运行到inf区域中，除非人们准确地设置了约束。我的问题:

1)我想念什么吗？是否有必要实现一些错误处理(对于erfc)？

2)您还有其他建议来加快代码或其他方法的速度吗？研究并行化for循环是否值得？

require(Rcpp)
require(RcppArmadillo)
require(microbenchmark)

#Example 1 : standard R::pnorm
src1 <- '
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const     arma::vec& sigma, int lt, int lg) {
int n = x.size();
arma::vec res(n);
for (int i=0; i<n; i++) {
   res(i) = R::pnorm(x(i),mu(i),sigma(i),lt,lg);
}
return wrap(res);
}
'

#Example 2: approximation with ercf
src2 <- '
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const    arma::vec& sigma, int lt, int lg) {
int n = x.size();
arma::vec res(n);
for (int i=0; i<n; i++) {
res(i) = 0.5 * erfc(-(x(i) - mu(i))/sigma(i) * M_SQRT1_2);
}
if (lt==0 & lg==0) {
   return wrap(1 - res);
}
if (lt==1 & lg==0) {
   return wrap(res);
}
if (lt==0 & lg==1) {
   return wrap(log(1 - res));
}
if (lt==1 & lg==1) {
   return wrap(log(res));
}
}
'

#some random numbers
xex  = rnorm(100,5,4)
muex = rnorm(100,3,1)
siex = rnorm(100,0.8,0.3)

#compile c++ functions
func1 = cppFunction(depends = "RcppArmadillo",code=src1) #R::pnorm
func2 = cppFunction(depends = "RcppArmadillo",code=src2) #ercf

#run with exemplaric data
res1 = func1(xex,muex,siex,1,0)
res2 = func2(xex,muex,siex,1,0)

# sum of squared errors
sum((res1 - res2)^2,na.rm=T)
# 6.474419e-32 ... very small

#benchmarking
 microbenchmark(func1(xex,muex,siex,1,0),func2(xex,muex,siex,1,0),times=10000)
#Unit: microseconds
#expr    min      lq     mean median     uq     max neval
#func1(xex, muex, siex, 1, 0) 11.225 11.9725 13.72518 12.460 13.617 103.654 10000
#func2(xex, muex, siex, 1, 0)  8.360  9.1410 10.62114  9.669 10.769 205.784 10000
#my machine: Ubuntu 14.04 LTS, i7 2640M 2.8 Ghz x 4, 8GB memory, RRO 3.2.0 based on version R 3.2.0

最佳答案

1)好吧，您确实应该使用R的pnorm()作为第0个示例。
不需要，您可以使用Rcpp接口(interface)。 R的pnorm()已经在R内部进行了很好的矢量化处理(即在C级别上)，因此可能比Rcpp甚至更快。它还具有覆盖NA，NaN，Inf等的情况的优势。

2)如果您正在谈论MLE，并且担心速度和准确性，则几乎可以肯定地应该使用对数，而不是pnorm()而不是dnorm()？

关于c++ - 计算 vector 正态分布的cdf的最快方法-R::pnorm vs erfc vs？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/30330569/