更新:根据其他人在源中找到的内容,我对此有误 - sum() 不排序。我在下面发现的一致性模式源于这样一个事实,即排序(如下面某些情况下所做的那样)和使用扩展精度中间值(如 sum() 所做的那样)会对精度产生类似的影响......
@user2357112 cmets 下面:
src/main/summary.c ... 不进行任何排序。 (添加到求和操作中会花费很多。)它甚至没有使用成对或补偿求和;它只是天真地将所有内容从左到右添加到 LDOUBLE 中(long double 或 double,取决于 HAVE_LONG_DOUBLE)。
我已经筋疲力尽地在 R 源代码中寻找这个(没有成功 - sum 很难搜索),但是 我可以通过实验证明在执行 sum() 时,R 对输入向量进行排序从最小到最大以最大限度地提高准确性;以下sum() 和Reduce() 结果之间的差异是由于使用了扩展精度。我不知道accu 做了什么……
set.seed(101)
vec <- runif(100, 0, 0.00001)
options(digits=20)
(s1 <- sum(vec))
## [1] 0.00052502325481269514554
使用Reduce("+",...) 只是按顺序添加元素。
(s2 <- Reduce("+",sort(vec)))
## [1] 0.00052502325481269514554
(s3 <- Reduce("+",vec))
## [1] 0.00052502325481269503712
identical(s1,s2) ## TRUE
?sum() 也说
在可能的情况下使用扩展精度累加器,但这取决于平台。
在 RcppArmadillo 中对已排序的向量执行此操作会得到与 R 中相同的答案;以原始顺序在向量上执行此操作会给出不同的答案(我不知道为什么;我的猜测是前面提到的扩展精度累加器,当数据未排序时,它会更多地影响数值结果)。
suppressMessages(require(inline))
code <- '
arma::vec ax = Rcpp::as<arma::vec>(x);
return Rcpp::wrap(arma::accu(ax));
'
## create the compiled function
armasum <- cxxfunction(signature(x="numeric"),
code,plugin="RcppArmadillo")
(s4 <- armasum(vec))
## [1] 0.00052502325481269525396
(s5 <- armasum(sort(vec)))
## [1] 0.00052502325481269514554
identical(s1,s5) ## TRUE
但正如 cmets 中所指出的,这并不适用于所有种子:在这种情况下,Reduce() 结果更接近sum()
的结果
set.seed(123)
vec2 <- runif(50000,0,0.000001)
s4 <- sum(vec2); s5 <- Reduce("+",sort(vec2))
s6 <- Reduce("+",vec2); s7 <- armasum(sort(vec2))
rbind(s4,s5,s6,s7)
## [,1]
## s4 0.024869900535651481843
## s5 0.024869900535651658785
## s6 0.024869900535651523477
## s7 0.024869900535651343065
我被难住了。我预计至少 s6 和 s7 是相同的......
我会指出,一般来说,当您的算法依赖于这些微小的数值差异时,您可能会感到非常沮丧,因为结果可能会因许多微小且可能超出的情况而有所不同-您的控制因素,例如您使用的特定操作系统、编译器等。