20

我很想知道 R 的平均函数使用什么算法。该算法的数值属性是否有一些参考？

我在 summary.c:do_summary() 中找到了以下 C 代码：

case REALSXP:
PROTECT(ans = allocVector(REALSXP, 1));
for (i = 0; i < n; i++) s += REAL(x)[i];
s /= n;
if(R_FINITE((double)s)) {
    for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
    s += t/n;
}
REAL(ans)[0] = s;
break;

这似乎是直截了当的意思：

for (i = 0; i < n; i++) s += REAL(x)[i];
s /= n;

然后它添加了我假设的数字校正，这似乎是与数据平均值的平均差异：

for (i = 0; i < n; i++) t += (REAL(x)[i] - s);
s += t/n;

我无法在任何地方追踪这个算法（平均值不是一个很好的搜索词）。

任何帮助将非常感激。

4

2 回答 2

15

我不确定这是什么算法，但 Martin Maechler 提到了West, 1979响应PR#1228的更新方法，该方法由 Brian Ripley 在 R-2.3.0 中实现。我在列出所使用的实际算法的源代码或版本控制日志中找不到参考。它cov.c在修订版 37389 和summary.c修订版 37393 中实施。

于 2013-07-25T19:52:20.590 回答
11

我相信 R 算法的工作原理如下。

平均值的第一个标准计算实际上是对代数平均值的估计，这是由于浮点误差（总和越远离被累积的元素，它会变得更糟）。

第二遍将元素与估计平均值的差异相加。应该没有净差​​，因为平均值两侧的值应该平衡，但我们有浮点误差。与平均值的差异仍然有可能出错，但这些应该小于元素与累积和之间的最差潜在差异（至少估计的平均值在值范围内的某个地方，而总和可能会逃脱它） . 除以 N 为您提供与平均值的平均差，然后您可以使用它来推动您的初始估计更接近真实平均值。您可以重复此操作以越来越接近，但在某些时候，计算与平均值的平均差异时的浮点误差会打败您。我猜一关就够近了。

这是我妻子向我解释的。

我不确定算法的来源是什么，我不确定这与其他方法相比如何，例如 Kahan 求和。我想我得做一些测试。

于 2013-07-26T14:14:49.090 回答