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

Question

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

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

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;

我无法在任何地方跟踪此算法(意味着不是一个很棒的搜索词).

任何帮助将非常感激.

Answer 1

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

Answer 2

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

平均值的第一个标准计算实际上是由于浮点误差导致的代数均值的估计(随着总和从积累的元素越远而变得越来越差).

第二遍将元素与估计平均值的差异相加.应该没有净差​​异,因为均值两边的值应该平衡,但我们有浮点误差.与均值的差异仍然存在误差的可能性,但是这些应该小于元素和累积和之间的最差电位差(至少估计的平均值存在于值范围内的某个位置,而求和可能会逃避它) .除以N可得出与平均值的平均差值,然后用它来推动初始估计值接近真实均值.您可以重复此操作以越来越近,但在某些时候,计算与平均值的平均差异时的浮点误差将使您失败.我猜一次传球足够接近.

这是我妻子向我解释的.

我不确定算法的来源是什么,我不确定这与其他方法相比如Kahan求和.我想我必须做一些测试.