计算向量上正态分布的 cdf 的最快方法 - R::pnorm vs erfc vs?

Posted 2023-02-16

【中文标题】计算向量上正态分布的 cdf 的最快方法 - R::pnorm vs erfc vs?

【英文标题】：Fastest way to compute the cdf of the Normal distribution over vectors - R::pnorm vs erfc vs?

【发布时间】：2015-05-19 15:59:27

【问题描述】：

我希望我重新措辞的问题现在符合 *** 的标准。请考虑以下示例。我正在编写一个对数似然函数，其中通过向量计算 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

【问题讨论】：

你读过this Rcpp Gallery post on subsetting吗？

如果你阅读了这篇文章，你为什么将犰狳函数应用于 Rcpp 类型？此外，您的示例既不完整也不可重复——而我们在 Rcpp Gallery 的文档既不完整也不可重复。我会从那里开始并进行相应的修改。

我们经常使用术语最小可重现示例。您的 codedump 既不是最小的，也不是（由于缺少输入数据）可重现的。

所以我完全重写了我的问题，现在，没有任何评论，它被搁置了吗？我猜之前投票删除它的人甚至没有审查新版本。

多个 cmet 教我如何提出我的问题，但没有对问题本身发表任何评论。伤心。

【参考方案1】：

1) 好吧，你真的应该使用 R 的 pnorm() 作为你的第 0 个例子。 你没有，你使用它的 Rcpp 接口。 R 的pnorm() 已经很好地在 R 内部进行了矢量化（即在 C 级别上），因此很可能是可比较的，甚至比 Rcpp 更快。此外，它确实具有覆盖 NA、NaN、Inf 等情况的优势。

2) 如果您谈论的是 MLE，并且您关心速度和准确性，那么您几乎可以肯定应该使用对数，也许不是pnorm()，而是dnorm()？

【讨论】：

感谢您的回复：1) 根据上面 Dirk Eddelbuettel 的评论，我假设 pnorm 的矢量化版本和我实际写出的代码是等价的。 2）我当然使用日志。也许我的语言很草率，应该写到这将成为对数似然函数的一部分。但是您能否概述一下我如何使用 dnorm 以及这会有什么优势？ MLE 是一个分层有序概率模型，在我阅读您的评论之前，我一直认为我需要 c.d.f.因此 pnorm ....

如果您确实需要概率模型，您确实需要pnorm()。 OTOH，许多统计学家确实更喜欢“默认” logit 而不是 probit ...一个原因是可能性变得更容易，包括它的梯度等。因此，对于您的问题，您可以先计算 logit 解决方案（它有一个“便宜的”对数可能性），然后将其解决方案用作更昂贵的概率的“初始起点”。我建议保留基于 pnorm() 的计算。速度损失小于 50%，那么您为什么要关心呢？我想说的是可靠性和可移植性更重要。