R使用什么算法来计算平均值？

Question

我很想知道R语言中mean函数使用的算法是什么。这个算法的数值特性有哪些参考资料吗？

我在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;

我无法在任何地方找到这个算法（mean不是一个很好的搜索词）。

任何帮助都将不胜感激。

- Zach

---

这只是一个旁白，但是'mean.R'如何调用'summary.c'呢？我不明白'.Internal(mean(x))'如何调用'summary.c'。感谢任何关于这两个文件如何连接的指针。如果这距离您的问题太远，对不起。我只是希望学习。 - Mark Miller

@MarkMiller：所有的.Internal调用都映射在src/main/names.c中。在该文件中搜索“mean”，您将看到调用它的C函数。然后，您可以搜索该C函数的源文件。请参见此问题。 - Joshua Ulrich

将此问题链接到r-devel：https://stat.ethz.ch/pipermail/r-devel/2013-July/067053.html - Brian Diggs

Answer 1

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

- Joshua Ulrich

---

谢谢你指引我正确的方向，我得去找一份这篇论文的副本。 - Zach

我认为那不是West的方法。我刚刚下载了这篇论文，West提出了一种单遍计算方差的方法 - R使用双遍方法计算平均值。 - hadley

我们在谈论这个修订版，对吧？链接 - hadley

@hadley：是的，那就是修订版。感谢您检查West算法。 - Joshua Ulrich

Answer 2

我认为 R 算法的工作原理如下。

第一个标准计算平均数实际上是代数平均数的估计，由于浮点误差（离累加元素越远，误差就会越大）。

第二个步骤对元素与估计平均值之差进行求和。由于平均值两侧的值应该相等，因此不应该存在净差异，但是我们有浮点误差。与平均值的差异仍然具有潜在的误差可能性，但这些误差应该小于元素与累加和之间最糟糕的潜在差异（至少估计的平均值位于值范围内，而总和可能会超出其范围）。除以 N 得到平均值与差距，然后使用它将初始估计值推近真实平均值。您可以重复此操作以逐渐接近真实平均值，但是在某个时候，计算平均值与差异的浮点误差将击败您。我想一次通行足够接近了。

这是我妻子告诉我的。

我不确定算法的来源在哪里，也不确定它与其他方法（如 Kahan 求和）相比如何。我想我需要做一些测试。

- Zach